library(survival)
library(stargazer)

d10merge<- read.csv("DavisWilf_JoiningTheClub_ReplicationData.csv")

	#Base sample, time to apply
		aa2 <- subset(d10merge, (ApplyYrI_NotYetApplied_a==1|ApplyYrI_Ind_a==1) 
				&end0>=1 & year>1947	& cty_ind==1); nrow(aa2) #3350
		length(unique(aa2$ccode))	#205 countries eligible to apply
		sum(aa2$ApplyYrI_Ind_a)		#185 countries apply
		nrow(subset(aa2, year==endyr_DWanalysis 
			& ApplyYrI_Ind_a==0))	# 20 countries never apply
		min(aa2$year); max(aa2$year)	#1948-2015 (though, base data through 2014 only!)

#####################################################################
#######
# Figure 1, left graph

	ctys <- subset(d10merge, styr_DWanalysis==year & (ApplyYrI_NotYetApplied_a==1 | ApplyYrI_Ind_a==1))
		nrow(ctys) #205
		head(ctys)

	timetoapply <- tta <- rep(NA, nrow(ctys))
	for(i in 1:nrow(ctys)){
		sub <- ctys[i,]
		timetoapply[i] <- ifelse(sub$ApplyYrI==0, sub$endyr_DWanalysis - sub$styr_DWanalysis + 1,
						sub$ApplyYrI - sub$styr_DWanalysis + 1)}
	sort(timetoapply)
		vals <- sort(unique(timetoapply))

	x <- 1:70
	count <- rep(0, length(x))
	h <- as.data.frame(cbind(x,count))
		head(h)
	for(i in 1:length(vals)){
		h$count[h$x==vals[i]] <- sum(timetoapply==vals[i])}

	B <- h$count

#	pdf("JOP_Figure1_left.pdf")
		barplot(B, col="grey")
			mids <- barplot(B, xlab="",ylab="",axes=F)
			axis(1, at=c(mids[1],mids[10],mids[20],mids[30],mids[40],mids[50],mids[60],mids[70]),
				labels=c("1","10","20","30","40","50","60","70"))
			axis(2, las=2)

		mtext("Years Between Eligibility and Application", side=1, line=2.5)
		mtext("Number of Countries", side=2, line=2.5)
		mtext("Distribution of Time to Apply", side=3, line=1,font=2,cex=1.5)
#	dev.off()
	
	##########################
	# "While almost 25 percent of all countries apply 
	# the first year they are eligible and 45 percent apply within 
	# the first five years of eligibility, others substantially 
	# delay application."

	#24.2% of countries apply in eligibility year 1
			sum(h$count)		#206 total countries
			sum(h$count[h$x==1])	#50 countries apply in year 1
		50/206 				# == 24.2%

	#45.6% of countries apply in eligibility year <=5
			sum(h$count)		#206 total countries
			sum(h$count[h$x<=5])	# 94 countries apply in year 1
		94/206 				# == 45.6%

#######
# Figure 1, right graph

	ctys <- 	subset(d10merge, NegotSample==1)
		nrow(ctys)	#101

	gatt <- subset(ctys, AccessionType=="GATT, Art33")
		gatt$NegSample_YrsToJoin
	wto <- subset(ctys, AccessionType!="GATT, Art33")
		wto$NegSample_YrsToJoin
	
		g <- length(gatt$NegSample_YrsToJoin)
		w <- length(wto$NegSample_YrsToJoin)
		g+w #101

	negottime <- nt <- rep(NA, nrow(ctys))
	for(i in 1:nrow(ctys)){
		sub <- ctys[i,]
		negottime[i] <- ifelse(sub$JoinYear==0, sub$endyr_DWanalysis-sub$ApplyYrI + 1,
						sub$JoinYear - sub$ApplyYrI + 1)}
	sort(negottime)	#1 - 32
		vals <- sort(unique(negottime))

	#GATT VALS
		ntGATT <- rep(NA, g)
		for(i in 1:g){
			sub <- gatt[i,]
			ntGATT[i] <- ifelse(sub$JoinYear==0, sub$endyr_DWanalysis-sub$ApplyYrI + 1,
						sub$JoinYear - sub$ApplyYrI + 1)}
		sort(ntGATT)	#1 - 32
			gvals <- sort(unique(ntGATT))

	#WTO VALS
		ntWTO <- rep(NA, w)
		for(i in 1:w){
			sub <- wto[i,]
			ntWTO[i] <- ifelse(sub$JoinYear==0, sub$endyr_DWanalysis-sub$ApplyYrI + 1,
						sub$JoinYear - sub$ApplyYrI + 1)}
		sort(ntWTO)	#1 - 29
			wvals <- sort(unique(ntWTO))

	x <- 1:35
	ctGATT <- rep(0, length(x))
	ctWTO <- rep(0, length(x))
	h <- as.data.frame(cbind(x,ctGATT,ctWTO))
		head(h)
	for(i in 1:length(gvals)){h$ctGATT[h$x==gvals[i]] <- sum(ntGATT==gvals[i])}
	for(i in 1:length(wvals)){h$ctWTO[h$x==wvals[i]] <- sum(ntWTO==wvals[i])}

	TOT <- h$TOT <- h$ctGATT+h$ctWTO
	G <- h$ctGATT
	W <- h$ctWTO

	#pdf("JOP_Figure1_right.pdf")
		barplot(B,col=c("black","white"))

	mids <- barplot(TOT,col="white", axes=F, ylab="",ylim=c(0,13),
		legend.text = c("GATT accession", "WTO accession"),
	        args.legend = list(x = "topright", fill=c("black","white"))
		)
		barplot(G, col=1, add=T,axes=F)
			axis(1, at=c(mids[1],mids[10],mids[20],mids[30],mids[35]),
				labels=c("1","10","20","30","35"))
			axis(2, las=2)
		mtext("Years From Application to Accession", side=1, line=2.5)
		mtext("Number of Countries", side=2, line=2.5)
		mtext("Distribution of Negotiation Time", side=3, line=1,font=2,cex=1.5)
	#dev.off()

#####################################################################
# Table 1

	nrow(simple <- subset(aa2, begin0==0))	#205
	simple$Art26_Strata_Dummy <- ifelse(simple$Art26_Strata==1,1,0)
	simple$YrsToApply <- ifelse(simple$ApplyYrI==0, simple$DWendyear - simple$year + 1, simple$ApplyYrI - simple$year + 1)
	simple$ApplyCensor <- ifelse(simple$ApplyYrI==0,0,1)
		sort(unique(simple$YrsToApply))
		
	#Formal Test of difference in survival times among groups of obs in StrataNegRig (1=founders, Art26; 2=Art33,WTO12)
		#	Testing the null hypothesis: ApplyTimes for groups in StrataCat2 are indistinguishable.
		survdiff(Surv(YrsToApply,ApplyCensor)~Art26_Strata_Dummy, rho=0,data=simple)	#default, rho=0... equivalent of the logrank test! S(t)^rho.
			#p=.00805 --> reject the null hypothesis; Strata groups have statistically different failure times.

	####DATA USED:. . . .		
	aa2 <- subset(d10merge, (ApplyYrI_NotYetApplied_a==1|ApplyYrI_Ind_a==1) 
				&end0>=1 & year>1947	& cty_ind==1); nrow(aa2) #3350

	#MODEL 1
	m1<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2
		+lnopentotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(aa2$polity2,aa2$lnopentotperc,aa2$lnnew_gdp_67,aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 2167
		length(unique(consol.1a.sample.in$ccode))	#countries = 144
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 137
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
		#CORRELATION MATRIX
			s1 <- consol.1a.sample.in
			hold1 <-data.frame(s1[,c("polity2", 
					"lnopentotperc",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd")])
			round(cor(hold1),3)

	#Model 2
	m2 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
			#Model 2, with time interaction for s3un100 <<RESULTS IN APPENDIX TABLE A.3, "Include Time Interaction">>
			m2_alt <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
				+polity2 +s3un100 +s3un100:log(end0) +allies_atop_ct_mbrs_zeroed
				+ lnopentotperc +gatttotperc
				+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
				+cluster(ccode)+strata(Art26_Strata),
				robust=TRUE,data=aa2); cox.zph(control.1a)
			summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1751
		length(unique(consol.1a.sample.in$ccode))	#countries = 134
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 121
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
		#CORRELATION MATRIX
			s2 <- consol.1a.sample.in
			hold2 <-data.frame(s2[,c("polity2", "s3un100", "allies_atop_ct_mbrs_zeroed",
					"lnopentotperc","gatttotperc",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd")])
			round(cor(hold2),3)

	#Model 3
	m3 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc + PercWorldMbrs
		+ trround_ind +trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1750
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
		#CORRELATION MATRIX
			s3 <- consol.1a.sample.in
			hold3 <-data.frame(s3[,c("polity2", "s3un100", "allies_atop_ct_mbrs_zeroed",
					"lnopentotperc","gatttotperc",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd",
					"Democratization_MP_ind_Polity",
					"colony_ind_gowa", "ptatrade_perc",
					"PercWorldMbrs", "trround_ind")])
			round(cor(hold3),3)

	#Model 4
	m4 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2_2yrchange +s3un100_2yrchange +allies_atop_ct_mbrs_2yrchange_zeroed
		+ opentotperc_2yrchange	+gatttotperc_2yrchange
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+ colony_ind_gowa +ptatrade_perc + PercWorldMbrs
		+ trround_ind + trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2_2yrchange, aa2$s3un100_2yrchange, 
					aa2$allies_atop_ct_mbrs_2yrchange_zeroed,
				aa2$opentotperc_2yrchange, aa2$gatttotperc_2yrchange,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1450
		length(unique(consol.1a.sample.in$ccode))	#countries = 91
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures =  78
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1950-2014
		#CORRELATION MATRIX
			s4 <- consol.1a.sample.in
			hold4 <-data.frame(s4[,c("polity2_2yrchange", "s3un100_2yrchange", 
					"allies_atop_ct_mbrs_2yrchange_zeroed",
					"opentotperc_2yrchange","gatttotperc_2yrchange",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd",
					"Democratization_MP_ind_Polity",
					"colony_ind_gowa", "ptatrade_perc",
					"PercWorldMbrs", "trround_ind")])
			round(cor(hold4),3)

	#Model 5
	negsampTV <- subset(d10merge, NegotSampleTV==1); 
	negsamp   <- subset(d10merge, NegotSampleTV==1 & ApplyYrI_Ind_a==1 ); 
			nrow(negsampTV)	#964
			nrow(negsamp)	#101	#101 total countries that apply under GATT Art 33 or WTO Art 12

	m5 <- control.1a<- coxph(Surv(NegSample_TV_begin0,NegSample_TV_end0,
						NegSample_TV_Censor_ieJoinInd) ~ 
		+ polity2 + s3un100 + allies_atop_ct_mbrs_zeroed
		+ lnopentotperc + gatttotperc
		+ lnnew_gdp_67 + lnpcgdp + ColdWarPd + ColdWarPd:log(NegSample_TV_end0)
		+cluster(ccode),robust=TRUE,data=negsampTV); cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				negsampTV$polity2, negsampTV$s3un100, negsampTV$allies_atop_ct_mbrs_zeroed,
				negsampTV$lnopentotperc, negsampTV$gatttotperc,
				negsampTV$lnnew_gdp_67,	negsampTV$lnpcgdp,
				negsampTV$ColdWarPd)
			sample<-cbind(negsampTV,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 680
		length(unique(consol.1a.sample.in$ccode))	#countries = 80
		sum(consol.1a.sample.in$NegSample_TV_Censor_ieJoinInd)	
									#failures =  62
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
		#CORRELATION MATRIX
			s5 <- consol.1a.sample.in
			hold5 <-data.frame(s5[,c("polity2", "s3un100", "allies_atop_ct_mbrs_zeroed",
					"lnopentotperc","gatttotperc",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd")])
			round(cor(hold5),3)

##########################################################################

#######
# Figure 2 -- left graph
#Estimated Probability of Application 

	#Model 3
	m3 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc + PercWorldMbrs
		+ trround_ind #+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); 
		cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1750
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
		#CORRELATION MATRIX
			s3 <- consol.1a.sample.in
			hold3 <-data.frame(s3[,c("polity2", "s3un100", "allies_atop_ct_mbrs_zeroed",
					"lnopentotperc","gatttotperc",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd",
					"Democratization_MP_ind_Polity",
					"colony_ind_gowa", "ptatrade_perc",
					"PercWorldMbrs", "trround_ind")])
			round(cor(hold3),3)

	#1st and 3rd quartile affinity values
		#SORTED IN-SAMPLE VALUES OF UN VOTING SIMILARITY
		hold <- s3[,c("statenme","year","s3un100")]
		hold[order(hold$s3un100),]

		# 1st and 3rd quartile values of in-sample affinity
			affinity[2]	#1st quartile affinity == -44.05
			affinity[5]	#3rd quartile affinity == +25.86
		ctys <- unique(hold$statenme)
		avgs3un100 <- rep(NA,length(ctys))
		for(i in 1:length(ctys)){
			sub <- subset(hold, statenme==as.character(ctys[i]))	
			avgs3un100[i] <- round(as.numeric(mean(sub$s3un100)),2)
			}
		h <- as.data.frame(cbind(as.character(ctys),avgs3un100))
			subset(h,avgs3un100== -43.97)	#Bulgaria
			subset(h,avgs3un100== 25.50)	#Bolivia

		cty <-subset(hold,statenme=="Bulgaria" & year>1956) #1981-86
			sort(cty$year)	#1948-1979
			mean(cty$s3un100)	#-47.56
		cty <-subset(hold,statenme=="Bolivia")	#1948-87
			sort(cty$year)	#1948-1979
			mean(cty$s3un100)	#25.499

	datapull<-as.data.frame(consol.1a.sample.in); nrow(datapull)
	affinity	<-as.numeric(summary(datapull$s3un100))
	ally		<-as.numeric(summary(datapull$allies_atop_ct_mbrs_zeroed))
	gatttr	<-as.numeric(summary(datapull$gatttotperc))
	cw		<-as.numeric(summary(datapull$ColdWarPd ))
	gdp		<-as.numeric(summary(datapull$lnnew_gdp_67))
	gdppc		<-as.numeric(summary(datapull$lnpcgdp))
	polity	<-as.numeric(summary(datapull$polity2))
	open		<-as.numeric(summary(datapull$lnopentotperc))
	colony	<-as.numeric(summary(datapull$colony_ind_gowa))
	round		<-as.numeric(summary(datapull$trround_ind ))
	mbrsperc	<-as.numeric(summary(datapull$PercWorldMbrs))
	ptadep	<-as.numeric(summary(datapull$ptatrade_perc))
	democratiz	<-as.numeric(summary(datapull$Democratization_MP_ind_Polity))

	affinity.surv<-data.frame(
		s3un100=c(affinity[2],affinity[5]),	#move from q1 to q3
		allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=rep(cw[4],2),
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=rep(polity[4],2),
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)


	par(mfrow=c(1,1))
		#time periods
	summary(survfit(control.1a, newdata=affinity.surv)[2])
		#data
	g <- summary(survfit(control.1a, newdata=affinity.surv)[2])
		g[6]	#2 columns - survival for values 1 and 2
			g[6][[1]]
			class(g[6])	#LIST
		h <- g[6][[1]]	
		h[,1]	#<- survival rates for first value
		h[,2]	#<- survival rates for first value

		1-h[,1]	#probability of failure, value 1
		1-h[,2]	#probability of failure, value 2
	#plus, add year 0 and end year for plotting(!)
	(y 	<- c(0,g[2][[1]]))	
	#plus, add beginning value, and additional final value(!)
	(survprobq1<- c( 1,   h[,1]))
	(failprobq1<- c( 0, 1-h[,1]))
	(survprobq3<- c( 1,   h[,2]))
	(failprobq3<- c( 0, 1-h[,2]))	

#	pdf("JOP_Figure2_left.pdf")
	par(mfrow=c(1,1),mar=(c(5, 6, 4, 2) + 0.1))	#default mar=c(5, 4, 4, 2) + 0.1
	plot(1:10,1:10, type="n",axes=F, 
		xlim=c(0,56), ylim=c(0,1), xlab="", ylab="")
	mtext("Estimated Probability of Applying \n UN Voting Similarity", side=3, line=1,font=2,cex=1.5)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Estimated Cumulative Probability\nof Having Applied", side=2, line=3)

	cbind(y,failprobq1,failprobq3)
		#plot horizontal lines for low values
		for(i in 1:(length(y)-1)){
			x1 <- y[i]
			x2 <- y[i+1]
			y1 <- failprobq1[i]
			y2 <- failprobq1[i]
			lines(c(x1,x2),c(y1,y2), lty=3, lwd=2)}
		#plot vertical lines for low values
		for(i in 1:(length(y)-1)){
			x1 <- y[i+1]
			x2 <- y[i+1]
			y1 <- failprobq1[i]
			y2 <- failprobq1[i+1]
			lines(c(x1,x2),c(y1,y2), lty=3, lwd=2)}
	cbind(y,failprobq1,failprobq3)
		#plot horizontal lines for high values
		for(i in 1:(length(y)-1)){
			x1 <- y[i]
			x2 <- y[i+1]
			y1 <- failprobq3[i]
			y2 <- failprobq3[i]
			lines(c(x1,x2),c(y1,y2), lty=1, lwd=2)}
		#plot vertical lines for high values
		for(i in 1:(length(y)-1)){
			x1 <- y[i+1]
			x2 <- y[i+1]
			y1 <- failprobq3[i]
			y2 <- failprobq3[i+1]
			lines(c(x1,x2),c(y1,y2), lty=1, lwd=2)}
			axis(1,at=c(0, 10,20,30,40,50,56),
					c("0","10","20","30","40","50","56"))
			axis(2,las=2)
		text(48,.35,"Low \n(1st quartile) \nUN voting\nsimilarity",cex=1.5)
		text(13,.85,"High \n(3rd quartile)\nUN voting\nsimilarity",cex=1.5)
#	dev.off()

#######
# Figure 2, right graph
##TIME TO APPLY FIRST DIFFERENCE CURVE


	##IDENTIFY FIRST DIFFERENCE SAMPLE THAT THE SIMULATION 
	## SHOULD DRAW FROM
	#Model 3
	m3 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc + PercWorldMbrs
		+ trround_ind #+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); 
		cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1750
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
		#CORRELATION MATRIX
			s3 <- consol.1a.sample.in
	sort(unique(s3$end0))#1-56
			hold3 <-data.frame(s3[,c("polity2", "s3un100", "allies_atop_ct_mbrs_zeroed",
					"lnopentotperc","gatttotperc",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd",
					"Democratization_MP_ind_Polity",
					"colony_ind_gowa", "ptatrade_perc",
					"PercWorldMbrs", "trround_ind")])
			round(cor(hold3),3)
	datapull<-as.data.frame(consol.1a.sample.in); nrow(datapull)
	affinity	<-as.numeric(summary(datapull$s3un100))
	ally		<-as.numeric(summary(datapull$allies_atop_ct_mbrs_zeroed))
	gatttr	<-as.numeric(summary(datapull$gatttotperc))
	cw		<-as.numeric(summary(datapull$ColdWarPd ))
	gdp		<-as.numeric(summary(datapull$lnnew_gdp_67))
	gdppc		<-as.numeric(summary(datapull$lnpcgdp))
	polity	<-as.numeric(summary(datapull$polity2))
	open		<-as.numeric(summary(datapull$lnopentotperc))
	colony	<-as.numeric(summary(datapull$colony_ind_gowa))
	round		<-as.numeric(summary(datapull$trround_ind ))
	mbrsperc	<-as.numeric(summary(datapull$PercWorldMbrs))
	ptadep	<-as.numeric(summary(datapull$ptatrade_perc))
	democratiz	<-as.numeric(summary(datapull$Democratization_MP_ind_Polity))

	affinity.surv<-data.frame(
		s3un100=c(affinity[2],affinity[5]),	#move from q1 to q3
		allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=rep(cw[4],2),
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=rep(polity[4],2),
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)
		
	#RUN THE SIMULATION 
	store.affCOUNTRY <- matrix(rep(NA,68000),nrow=68, ncol=1000)
	for(k in 1:1000){
		datapull<-as.data.frame(consol.1a.sample.in)
		indx<-sample(1:nrow(datapull),5000,replace=TRUE)
		b.data<-consol.1a.sample.in[indx,]

	fitcty	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc + PercWorldMbrs
		+ trround_ind #+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=b.data)
	
		temp <- summary(survfit(fitcty, newdata=affinity.surv))
		d <- data.frame(first = temp$surv[, 1], second = temp$surv[,2], strata = temp$strata, time = temp$time)
		d.subset <- d[d$strata == levels(d$strata)[2], ]
		d.diff <- d.subset[,2] - d.subset[,1]
		store <- rep(NA, 68)
		for(i in 1:68){
			if (i %in% d.subset$time) {
			     store[i] <- d.diff[d.subset$time == i]
			  } else if (i == 1) {
			     store[1] <- 0
			  } else {
			     store[i] <- store[i-1]
			}
			}
		store.affCOUNTRY[,k]<-store
		}
	hold<-store.affCOUNTRY
	hold
	hold5<-hold 	
	#	OR, can load previous set of estimates
	#hold5 <- read.csv("FirstDifference,SimulationData,Jop,Model3.csv")
	head(hold5)

	#TAKE AVERAGES AND 2.5% AND 97.5% OBSERVATIONS
	ctyest.diff<-ctyci.high<-ctyci.lo<-rep(NA,68)
	for(i in 1:68){
		ctyest.diff[i]<-mean(as.numeric(hold5[i,2:1001]),na.rm=T)
		ctyci.high[i]<-as.numeric(quantile(hold5[i,2:1001],.975,na.rm=T))
		ctyci.lo[i]<-as.numeric(quantile(hold5[i,2:1001],.025,na.rm=T))
		}
	fd<-as.data.frame(cbind(ctyest.diff,ctyci.lo,ctyci.high))
		head(fd)		

	#PLOT THE FIRST DIFFERENCE CURVE
#pdf("JOP_Figure2_right.pdf")
			#par(mfrow=c(1,1))
	par(mfrow=c(1,1),mar=(c(5, 6, 4, 2) + 0.1))	#default mar=c(5, 4, 4, 2) + 0.1
	plot(1:56, -(fd$ctyest.diff)[1:56], ylim=c(-.01,.50),xlim=c(0,56),
		ylab="", xlab="",
		type="l", lwd=2,axes=F)
	axis(1,at=c(0, 10,20,30,40,50,56),
					c("0","10","20","30","40","50","56"))
	axis(2,las=2)
	abline(h=0)
	lines(1:68, abs(fd$ctyci.high), lty=2, lwd=1.5)
	lines(1:68, abs(fd$ctyci.lo), lty=2, lwd=1.5)
	mtext("First Difference \n UN Voting Similarity", side=3, line=1,font=2,cex=1.5)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Cumulative Probability of Having Applied\nEstimated Difference, 1st to 3rd quartile", side=2, line=3)
#dev.off()

	#MAX DIFERENCE IS 32.25 % in year 39 / 1986
		cbind(1:68, 1948:(1948+67))
	h <- cbind(EligibleYr=fd[,1],CalendarYr=1947+fd[,1],FirstDiff=abs(fd$ctyest.diff))
		subset(h, h[,3]==max(h[,3]))


##########################################################################

#######
# Appendix Table A.1
##Countries in sample

	#SAMPLE THAT ENTERS INTO BASELINE MODEL (TABLE 1, MODEL 1) ANALYSIS
		consol.1a.sample <- complete.cases(aa2$polity2,aa2$lnopentotperc,aa2$lnnew_gdp_67,aa2$lnpcgdp,aa2$ColdWarPd)
		sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
			nrow(consol.1a.sample.in)			#obs = 2167
			length(unique(consol.1a.sample.in$ccode))	#countries = 144
			sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 137
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
	#Table 1  Model 1 sample, All countries
		sort(unique(consol.1a.sample.in$statenme))
	#Table 1  Model 1 sample, GATT, Art 26 Eligible
		sort(unique(subset(consol.1a.sample.in, Art26_Strata==1)$statenme))
	#Table 1  Model 1 sample, Founders
		sort(unique(subset(consol.1a.sample.in, AppType_Founder==1)$statenme))

	#COUNTRIES IN NEGOTIATION ANAYSIS (TABLE 1, MODEL 5) SAMPLE
	negsampTV <- subset(d10merge, NegotSampleTV==1); 
			nrow(negsampTV)	#964
	consol.1a.sample <- complete.cases(
				negsampTV$polity2, negsampTV$s3un100, negsampTV$allies_atop_ct_mbrs_zeroed,
				negsampTV$lnopentotperc, negsampTV$gatttotperc,
				negsampTV$lnnew_gdp_67,	negsampTV$lnpcgdp,
				negsampTV$ColdWarPd)
			sample<-cbind(negsampTV,consol.1a.sample);Table1Model5Sample<-consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(Table1Model5Sample)			#obs = 680
		length(unique(Table1Model5Sample$ccode))	#countries = 80
		sum(Table1Model5Sample$NegSample_TV_Censor_ieJoinInd)	
									#failures =  62
			min(Table1Model5Sample$year)
			max(Table1Model5Sample$year) 	#years = 1948-2014
	#Table 1  Model 5 sample, All countries
		sort(unique(Table1Model5Sample$statenme))


#######
# Appendix Table A.2
##Reduced models that maximize in-sample country count

	s1<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+lnnew_gdp_67+lnpcgdp +ColdWarPd
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)	#n=3199
		#SAMPLE
			consol.1a.sample <- complete.cases(aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnpcgdp)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 3111
		length(unique(consol.1a.sample.in$ccode))	#countries = 194
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 180
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	s2<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+lnopentotperc
		+lnnew_gdp_67+lnpcgdp +ColdWarPd
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)	
		#SAMPLE
			consol.1a.sample <- complete.cases(aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnopentotperc,
							aa2$lnpcgdp)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 2467
		length(unique(consol.1a.sample.in$ccode))	#countries = 165
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 162
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	s3 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2
		+lnnew_gdp_67		 +lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)	
		#SAMPLE
			consol.1a.sample <- complete.cases(aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnpcgdp,
							aa2$polity2)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 2556
		length(unique(consol.1a.sample.in$ccode))	#countries = 158
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 148
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	s4 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2
		+lnopentotperc
		+lnnew_gdp_67	+lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)	
		#SAMPLE
			consol.1a.sample <- complete.cases(aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnpcgdp,aa2$lnopentotperc,
							aa2$polity2)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 2167
		length(unique(consol.1a.sample.in$ccode))	#countries = 144
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 137
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	s5<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+s3un100
		+s3un100:log(end0)
		+lnnew_gdp_67+lnnew_gdp_67:log(end0)	+lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)	
			#SAMPLE
			consol.1a.sample <- complete.cases(aa2$s3un100,
							aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnpcgdp)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 2627
		length(unique(consol.1a.sample.in$ccode))	#countries = 176
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 158
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	s6<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+s3un100
		+s3un100:log(end0)
		+lnopentotperc
		+lnnew_gdp_67	+lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)	
			#SAMPLE
			consol.1a.sample <- complete.cases(aa2$s3un100,
							aa2$lnopentotperc,
							aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnpcgdp)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 2185
		length(unique(consol.1a.sample.in$ccode))	#countries = 152
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 145
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	s7<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+allies_atop_ct_mbrs_zeroed
		+lnnew_gdp_67	+lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)
			#SAMPLE
			consol.1a.sample <- complete.cases(aa2$allies_atop_ct_mbrs_zeroed,
							aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnpcgdp)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 3111
		length(unique(consol.1a.sample.in$ccode))	#countries = 194
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 180
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014


	s8<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+allies_atop_ct_mbrs_zeroed
		+lnopentotperc
		+lnnew_gdp_67	+lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)	
			#SAMPLE
			consol.1a.sample <- complete.cases(aa2$allies_atop_ct_mbrs_zeroed,
							aa2$lnopentotperc,
							aa2$ColdWarPd,
							aa2$lnnew_gdp_67,
							aa2$lnpcgdp)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 2467
		length(unique(consol.1a.sample.in$ccode))	#countries = 165
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 162
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014



#######
# Appendix Table A.3
##Table 1, Model 2, Economic Variable Significance 
##	and Proportional Hazard Assumption (PHA)

	#"Table 1, Model 2"
	m2 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1751
		length(unique(consol.1a.sample.in$ccode))	#countries = 134
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 121
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
		#CORRELATION MATRIX
			s2 <- consol.1a.sample.in
			hold2 <-data.frame(s2[,c("polity2", "s3un100", "allies_atop_ct_mbrs_zeroed",
					"lnopentotperc","gatttotperc",
					"lnnew_gdp_67","lnpcgdp","ColdWarPd")])
			round(cor(hold2),3)

	#"Exclude geopolitical variables"
	#(Use sample from Table 1, Model 2)
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);T1M2<-consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
	m2_alt1 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=T1M2); cox.zph(control.1a)
	summary(control.1a)

	#"Include time interaction"
	#Model 2, with time interaction for s3un100 <<RESULTS IN APPENDIX TABLE A.3, "Include Time Interaction">>
	m2_alt2 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
				+polity2 +s3un100 +s3un100:log(end0) +allies_atop_ct_mbrs_zeroed
				+ lnopentotperc +gatttotperc
				+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
				+cluster(ccode)+strata(Art26_Strata),
				robust=TRUE,data=aa2); cox.zph(control.1a)
			summary(control.1a)

#######
# Number of accessions via Article 26.
	#	unique(d10$AccessionType)
	length(unique(subset(d10, AccessionType=="GATT, Art26:5")$ccode)) #64
	length(unique(subset(d10, AppType_Gatt26==1)$ccode)) #64

#######
# Number of potential jurisdictions / states in the system (fn.9)

	#205 potential applicants
		aa2 <- subset(d10merge, (ApplyYrI_NotYetApplied_a==1|ApplyYrI_Ind_a==1) 
			&end0>=1 & year>1947	& cty_ind==1)
	length(unique(aa2$ccode))	#205

#######
#	"Of 205 jurisdictions that could potentially apply, 185 do.

	#205 potential applicants
		aa2 <- subset(d10merge, (ApplyYrI_NotYetApplied_a==1|ApplyYrI_Ind_a==1) 
			&end0>=1 & year>1947	& cty_ind==1)
	length(unique(aa2$ccode))	#205

	#185 observed applicants
		applicants <- subset(aa2, ApplyYrI_Ind_a==1)		
	nrow(applicants)	  	#185

	#20 jurisdictions never apply
		sub <- subset(aa2, (year==endyr_DWanalysis & ApplyYrI_NotYetApplied_a==1))
		nrow(sub)	#20
	sub[,c("statenme","endyr_DWanalysis","year")]

#######
#	Footnote xx, Stratification by Article 26:5 is appropriate. . . 

	#Formal Test of difference in survival times among groups of obs in StrataNegRig (1=founders, Art26; 2=Art33,WTO12)
	#	Testing the null hypothesis: 
	#	ApplyTimes for groups in StrataCat2 are indistinguishable.

	sub <- subset(aa2, ApplyYrI_Ind_a==1 |(year==endyr_DWanalysis & ApplyYrI_NotYetApplied_a==1))
 		nrow(sub)	#205
		sum(sub$Art26_Strata==1)	#n= 76 ctys are Article 26:5 eligible
		sum(sub$Art26_Strata==2)	#n=129 ctys non-Article 26:5 eligible

		#Censor indicator and status
		sum(sub$ApplyYrI_Ind_a==1)	#185 applicants
		sum(sub$ApplyYrI_Ind_a==0)	# 20 censored countries

		#Identify Years to Apply
		sub$YrsToApply <- ifelse(sub$year==2015, 
			2014 - sub$styr_DWanalysis + 1,
			sub$year - sub$styr_DWanalysis + 1)

	survdiff(Surv(YrsToApply,ApplyYrI_Ind_a)~Art26_Strata, rho=-1,data=sub)	#rho<0 gives weight to the later part of survival curve.
		#p=7.98*10^-8 --> reject the null hypothesis; Stratified groups have different failure times.
	survdiff(Surv(YrsToApply,ApplyYrI_Ind_a)~Art26_Strata, data=sub)	#default, rho=0... equivalent of the logrank test! S(t)^rho.
		#p=.00074 --> reject the null hypothesis; Stratified groups have different failure times.
	survdiff(Surv(YrsToApply,ApplyYrI_Ind_a)~Art26_Strata, data=sub, rho=1) #rho>0 gives weight to the first part of survival curve.
		#p=.260 --> cannot reject the null hypothesis; Stratified groups 
		#				do not have distinct failure times in a statistically significant way.

#######
# Table A.4, 
# 	Robustness, Table 1 of 3
	
	#Table A.4,Model (1)
	# Cheibub Vreeland and Ghandi
	#	Democracy
	rob1 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+democracy
		+s3un100
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+ trround_ind +trround_ind:log(end0)
		+ PercWorldMbrs + PercWorldMbrs:log(end0)
		+ ptatrade_perc + colony_ind_gowa +Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.4,Model (2)
	# Freedom House Democracy
	rob2 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+fhval_flip_7dem
		+s3un100
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		#+ PercWorldMbrs
			#+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.4,Model (3)
	#Distance from mean member polity score
	rob3<- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+distfrom_mbr_mean_polity +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc + PercWorldMbrs
		+ trround_ind 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
			consol.1a.sample <- complete.cases(
				aa2$distfrom_mbr_mean_polity, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);samplein<-subset(sample, consol.1a.sample == TRUE)
			nrow(samplein)	#1750

	#Table A.4,Model (4)
	# Ideal point
	rob4 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 
		+ideal1
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+ trround_ind 
			#+trround_ind:log(end0)
		+ PercWorldMbrs
			+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.4,Model (5)
	#UNGA Important Votes
	rob5 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +allies_atop_ct_mbrs_zeroed
		+s3unimportant100
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc + ptatrade_perc:log(end0)#+ PercWorldMbrs
		+ trround_ind 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.4,Model (6)
	# Small Group (5-Country) Affinity Average
	rob6 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 
		+SmallGpAffinity_Avg100
		+allies_atop_ct_mbrs_zeroed
			#+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		#	+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.4,Model (7)
	#ATOP not extended
	# end sample in 2012 (ally end year)
	rob7 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		#+ lnopentotperc +gatttotperc
		+ lnopentotperc 
		#	+ lnopentotperc:log(end0)
		+gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		#+ColdWarPd:log(end0)
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		#	+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=subset(aa2, year<=2012
						)); cox.zph(control.1a)
	summary(control.1a)

	#Table A.4,Model (8)
	#Ally Weighted (1 of 3)
	#	weighted by trade + full model
	rob8 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	
			+MbrAllyTrade_Perc_MbrTrade
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc + PercWorldMbrs+ PercWorldMbrs:log(end0)
		+ trround_ind +trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.4,Model (9)
	#Ally Weighted (2 of 3)
	#	weighted by trade + reduced model
	rob9 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+MbrAllyTrade_Perc_MbrTrade
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)	
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$MbrAllyTrade_Perc_MbrTrade,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 3111
		length(unique(consol.1a.sample.in$ccode))	#countries = 194
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 180
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	#Table A.4,Model (10)
	#Ally Weighted (3 of 3)
	#	weighted by GDP + reduced model
	rob10 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+MbrAllyGdp_Perc_MbrGdp
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$MbrAllyGdp_Perc_MbrGdp,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 3111
		length(unique(consol.1a.sample.in$ccode))	#countries = 194
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 180
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

#######
# Table A.5, 
#	Robustness Checks 2 of 3
	
	#Table A.5,Model (1)
	#Article 26 subsample (1 of 2)
	# WITH s3un100, though it violated PHA
	rob1 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc 
		+ PercWorldMbrs
		+ trround_ind 
		+cluster(ccode),
		robust=TRUE,data=subset(aa2, year<1995 & Art26_Strata==1)); cox.zph(control.1a)
		summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnoptentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);
			consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE
					& year<1995 & Art26_Strata==1)
		nrow(consol.1a.sample.in)			#obs = 406
		length(unique(consol.1a.sample.in$ccode))	#countries = 46
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 39
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-1994

	#Table A.5, Model (2)
	#Article 26 subsample (2 of 2)
	# EXCLUDE s3un100 (with same substantive results)
	# use same sample as the analysis in Table A.5, Model (1)
	rob2 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc 
		+ PercWorldMbrs
		+ trround_ind 
		+cluster(ccode),
		robust=TRUE, data=consol.1a.sample.in); cox.zph(control.1a)
		summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnoptentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd)
			sample<-cbind(aa2,consol.1a.sample);
			consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE
					& year<1995 & Art26_Strata==1)
		nrow(consol.1a.sample.in)			#obs = 406
		length(unique(consol.1a.sample.in$ccode))	#countries = 46
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 39
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-1994


	#Table A.5,Model (3)
	# No stratification of Art 26 eligibility
	# + control for Article 26 eligibility
	aa2$Art26_Strata_Dummy <- ifelse(aa2$Art26_Strata==1,1,0)
	rob3 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+Art26_Strata_Dummy
		+cluster(ccode),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
	
	#Table A.5,Model (4)
	# BRESLOW TIES
	rob4 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2,method="breslow"); cox.zph(control.1a)
	summary(control.1a)

	#Table A.5,Model (5)
	# TIME TO JOIN (alternative dependent variable)
	ttj <- subset(d10merge, (JoinYear==0 | year<=JoinYear) & end0>=1 & year>1947	& cty_ind==1); 
			nrow(ttj) #4221
	rob5 <-control.1a<- coxph(Surv(begin0,end0,JoinYear_Ind) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 
			+lnnew_gdp_67:log(end0)
		+lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=ttj); cox.zph(control.1a)
	summary(control.1a)

	#Table A.5,Model (6)
	#Calendar Year perspective
	rob6 <-control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.5,Model (7)
	# GATT PERIOD ONLY
	rob7 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		#+ColdWarPd:log(end0)
		+ trround_ind 
		#	+trround_ind:log(end0)
		+ PercWorldMbrs
		#	+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=subset(aa2,year<=1994)); cox.zph(control.1a)
	summary(control.1a)

	#Table A.5,Model (8)
	# EXCLUDE FOUNDING MEMBERS
	rob8 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		#+ColdWarPd:log(end0)
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		#	+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=subset(aa2,AppType_Founder!=1)); cox.zph(control.1a)
	summary(control.1a)

	#Table A.5,Model (9)
	# EXCLUDE EARLY JOINERS
	rob9 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		#+ColdWarPd:log(end0)
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		#	+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=subset(aa2,JoinYear==0 | JoinYear>1953)); cox.zph(control.1a)
	summary(control.1a)

	#Table A.5,Model (10)
	# WTO period control
	rob10 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+wtopd
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		#+ColdWarPd:log(end0)
		+ trround_ind 
		#	+trround_ind:log(end0)
		+ PercWorldMbrs
			+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

#######
# Table A.6, 
# 	Robustness, Table 3 of 3
	
	#Table A.6,Model (1)
	# Cold War Intx (1 of 3)
	#Cold War Pd * Polity Score
	rob1 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		#+ColdWarPd:log(end0)
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		#	+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
			+ColdWarPd:polity2
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.6,Model (2)
	# Cold War Intx (2 of 3)
	#Cold War Pd * UN Voting Similarity
	rob2 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
			+ColdWarPd:s3un100
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.6,Model (3)
	# Cold War Intx (3 of 3)
	#Cold War Pd * Ally Member Count
	rob3 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		#+ColdWarPd:log(end0)
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		#	+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
			+ColdWarPd:allies_atop_ct_mbrs_zeroed
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.6, Model (4)
	# New IGO memberships
	rob4 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+IgoChange
		+ trround_ind 
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.6, Model (5)
	# Anti-Dumping (AD) legislation
	rob5 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		adlaw
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.6, Model (6)
	# Intx: GATT/WTO Member Trade * PTA Trade
	rob6 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
			+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
			+gatttotperc:ptatrade_perc
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.6, Model (7)
	# Intx: GATT/WTO Member Trade * Democracy
	rob7 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
			+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
			+gatttotperc:polity2
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)

	#Table A.6, Model (8)
	# Openness Decomposition
	#Export and Import openness
	rob8 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed	
		+ lnopenXperc +lnopenMperc
		+gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
			+ PercWorldMbrs:log(end0)
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
		summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopenXperc, aa2$lnopenMperc,
				aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1750
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	#Table A.6, Model (9)
	# GATT/WTO Trade Decomposition
	#Export and Import GATT/WTO Trade Perc
	rob9 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gattXperc + gattMperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
			+trround_ind:log(end0)
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2); cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, #aa2$gatttotperc,
				aa2$gattXperc,aa2$gattMperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1695
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

#######
# Table A.7, 
#  Cold War Period Subset (1948-1991)

	#Cold War Period Sample
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);
		cw1<-	consol.1a.sample.in<-subset(sample, year<=1991 & consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1497
		length(unique(consol.1a.sample.in$ccode))	#countries = 116
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 76
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-1991

	# Table A.7, Model (1)
	cw1_1 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

	# Table A.7, Model (2)
	cw1_2 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 #+s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

	# Table A.7, Model (3)
	cw1_3 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+s3un100
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

	# Table A.7, Model (4)
	cw1_4 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

	# Table A.7, Model (5)
	cw1_5 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

	# Table A.7, Model (6)
	cw1_6 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+polity2 
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

	# Table A.7, Model (7)
	cw1_7 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+s3un100
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

	# Table A.7, Model (8)
	cw1_8 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw1); cox.zph(control.1a)
	summary(control.1a)

#######
# Table A.8 
#	Post-Cold War Period Subset (1992+)

	#Post-Cold War Period Sample
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);
		cw0<-	consol.1a.sample.in<-subset(sample, year>1991 & consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 253
		length(unique(consol.1a.sample.in$ccode))	#countries = 52
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 44
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1992-2014

	# Table A.8, Model (1)
	cw0_1 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
		+ ptatrade_perc:log(end0)
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)

	# Table A.8, Model (2)
	cw0_2 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2 
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
				+ ptatrade_perc:log(end0)
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)

	# Table A.8, Model (3)
	cw0_3 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+s3un100
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
				+ ptatrade_perc:log(end0)
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)

	# Table A.8, Model (4)
	cw0_4 <- control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
			+gatttotperc:log(end0)
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ ptatrade_perc 
				+ ptatrade_perc:log(end0)
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)

	# Table A.8, Model (5)
	cw0_5 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)

	# Table A.8, Model (6)
	cw0_6 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+polity2 
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)

	# Table A.8, Model (7)
	cw0_7 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+s3un100
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)

	# Table A.8, Model (8)
	cw0_8 <- control.1a<- coxph(Surv(begin1948,end1948,ApplyYrI_Ind_a) ~ 
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp 
		+ColdWarPd 
		+ trround_ind 
		+ PercWorldMbrs
		+ ptatrade_perc 
		+ colony_ind_gowa
		+Democratization_MP_ind_Polity 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=cw0); cox.zph(control.1a)
	summary(control.1a)


#######
# Figure A.1
#	Hypothetical Trade from Earlier Joining
grt <- read.dta("GRT_IO_2007.csv")
	head(grt)

	###CREATE INDICATORS FOR 'nonppt', 'nonmbrppt', 'mbr'
			#gatt_1 = gatt ppt status of cty 1
			#gatt_2 = gatt ppt status of cty 2
		#PPT STATUSES:
			# FORMAL MEMBERS == orig, art33, art26:5, wto
			# NONMBR PPTS == col, df, prov
			# OUT == non-ppt
	grt$mbr_ind_1 <- ifelse(	(grt$gatt_1=="orig" |
							grt$gatt_1=="art33" |
							grt$gatt_1=="art26:5" |
							grt$gatt_1=="wto" ), 1, 0)	#unique(grt[,c("gatt_1","mbr_ind_1")])
	grt$nonmbrppt_ind_1 <- ifelse(	(grt$gatt_1=="col" |
							grt$gatt_1=="df" |
							grt$gatt_1=="prov" ), 1, 0)	#unique(grt[,c("gatt_1","nonmbrppt_ind_1")])
	grt$nonppt_ind_1 <- ifelse(	(grt$gatt_1=="out"), 1, 0)		#unique(grt[,c("gatt_1","nonppt_ind_1")])
	grt$mbr_ind_2 <- ifelse(	(grt$gatt_2=="orig" |
							grt$gatt_2=="art33" |
							grt$gatt_2=="art26:5" |
							grt$gatt_2=="wto" ), 1, 0)	#unique(grt[,c("gatt_2","mbr_ind_2")])
	grt$nonmbrppt_ind_2 <- ifelse(	(grt$gatt_2=="col" |
							grt$gatt_2=="df" |
							grt$gatt_2=="prov" ), 1, 0)	#unique(grt[,c("gatt_2","nonmbrppt_ind_2")])
	grt$nonppt_ind_2 <- ifelse(	(grt$gatt_2=="out"), 1, 0)		#unique(grt[,c("gatt_2","nonppt_ind_2")])


#######################
# Identify countries of interest. . . 
	unique(grt[,c("name_1","tomz_1")])	
		#tomz_1
		#Mexico, 411
		#China, 621

	unique(grt[,c("name_2","tomz_2")])	
		#tomz_1
		#Mexico, 411
		#China, 621

	x<-621	
	t<-subset(grt, tomz_1==x); unique(t[order(t$year),c("name_1","year","gatt_1")])	
	#China
	#	1950; 2001-03, 	mbr		('orig')
	#	1951--2000,		nonppt	('out')
	#	2001-03, 		mbr		('wto')

	x<-411	
	t<-subset(grt, tomz_1==x); unique(t[order(t$year),c("name_1","year","gatt_1")])	
	#Mexico
	#	1946--85, 	nonppt	('out')
	#	1986--2004,	mbr		('art 33' and 'wto')


###AGGREGATE TRADE. . . by dyad type . . .  FOR CTYS of interest

	#imports = natural log of M by cty 1 FROM cty 2, in constant 1967 USD
	grt$imports_nonlog <- exp(grt$imports)

###AGGREGATE TRADE
#CHINA
	(yrs <- sort(unique(grt$year))) #1946-2004

	ctyX <- subset(grt, tomz_2 == 621)
	(Xyrs <- as.data.frame(sort(unique(ctyX$year)))) 	#1950-2003
	names(Xyrs) <- "yr"
	Xyrs$trade_withMbrs <-rep(0, length(Xyrs))
	Xyrs$trade_withNonMbrPpts<-rep(0, length(Xyrs))
	Xyrs$trade_withNonPpts<-rep(0, length(Xyrs))
		head(Xyrs)
	for(i in 1:nrow(Xyrs)){
		sub <- subset(ctyX, year==Xyrs$yr[i])
		Xyrs$trade_withMbrs[i] <- ifelse(sum(sub$mbr_ind_1==1)==0,0,sum(subset(sub,mbr_ind_1==1)$imports_nonlog))
		Xyrs$trade_withNonMbrPpts[i] <- ifelse(sum(sub$nonmbrppt_ind_1==1)==0,0,sum(subset(sub,nonmbrppt_ind_1==1)$imports_nonlog))
		Xyrs$trade_withNonPpts[i] <- ifelse(sum(sub$nonppt_ind_1==1)==0,0,sum(subset(sub,nonppt_ind_1==1)$imports_nonlog))
		}
	head(Xyrs)
		M 	<- Xyrs[,c("yr","trade_withMbrs")]
		NMP 	<- Xyrs[,c("yr","trade_withNonMbrPpts")]
		NP	<- Xyrs[,c("yr","trade_withNonPpts")]

		M1 <- Xyrs
	M1$yrM1 <- M1$yr+1		
			head(M1)
		subM <- M1[,c("yrM1","trade_withMbrs")]; names(subM)[2]<-"twM_m1"; head(subM)
		subNMP <- M1[,c("yrM1","trade_withNonMbrPpts")]; names(subNMP)[2]<-"twNMP_m1"; head(subNMP)
		subNP <- M1[,c("yrM1","trade_withNonPpts")]; names(subNP)[2]<-"twNP_m1"; head(subNP)

	#Export data, trade with formal members
	(Xm <- merge(M, subM, by.x="yr", by.y="yrM1", all.x=T))
			Xm$tradeincrease <- Xm[,2] - Xm[,3]; Xm
			Xm$tradeincrease[1] <- 0; Xm

			Xm$adj <- rep(0,nrow(Xm))				#hyp / obs
		Xm$adj[Xm$yr<=1988] 			<- 1.22 / 1.22	#(NP-M) 	/ (NP-M)
		Xm$adj[Xm$yr>=1989 & Xm$yr<=2000] 	<- 1.41 / 1.22	#(**M**-M) 	/ (NP-M)
		Xm$adj[Xm$yr>=2001] 			<- 1.41 / 1.41	#(M-M)	/ (M-M)
			Xm$adjustment <- rep(0,nrow(Xm))
		Xm$adjustment[Xm$yr<=1988] 			<- "=1.22/1.22"
		Xm$adjustment[Xm$yr>=1989 & Xm$yr<=2000] 	<- "=1.41/1.22"
		Xm$adjustment[Xm$yr>=2001] 			<- "=1.41/1.41"

		Xm$HypTr <- rep(0,nrow(Xm))
			Xm$HypTr[1] <- Xm$trade_withMbrs[1]
		for(i in 2:nrow(Xm)){
				(year <- Xm$yr[i])
				(year_m1 <- ifelse(i==1,Xm$yr[1],Xm$yr[i]-1))
			(sub <- subset(Xm, yr==year))
			(sub_m1 <- subset(Xm, yr==year_m1))
		(Xm$HypTr[i] <- ifelse(sub$trade_withMbrs==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	dm1 <- Xm[,c("yr","trade_withMbrs","HypTr")];names(dm1)[3]<-"hyptr_xdata_M"; head(dm1)

	#Export data, trade with non-mbr ppts
	(Xnmp <- merge(NMP,subNMP, by.x="yr", by.y="yrM1", all.x=T))
			Xnmp$tradeincrease <- Xnmp[,2] - Xnmp[,3]; Xnmp
			Xnmp$tradeincrease[1] <- 0; Xnmp

			Xnmp$adj <- rep(0,nrow(Xnmp))					#hyp 		/ obs
		Xnmp$adj[Xnmp$yr<=1988] 			<- 1.19 / 1.19	#(NP-NMP) 	 / (NP-NMP)
		Xnmp$adj[Xnmp$yr>=1989 & Xnmp$yr<=2000] 	<- 1.46 / 1.19	#(**M**-NMP) / (NP-NMP)
		Xnmp$adj[Xnmp$yr>=2001] 			<- 1.46 / 1.46	#(M-NMP)	 / ( M-NMP)
			Xnmp$adjustment <- rep(0,nrow(Xnmp))
		Xnmp$adjustment[Xnmp$yr<=1988] 			<- "=1.19/1.19"
		Xnmp$adjustment[Xnmp$yr>=1989 & Xnmp$yr<=2000] 	<- "=1.46/1.19"
		Xnmp$adjustment[Xnmp$yr>=2001] 			<- "=1.46/1.46"

		Xnmp$HypTr <- rep(0,nrow(Xnmp))
			Xnmp$HypTr[1] <- Xnmp$trade_withNonMbrPpts[1]
		for(i in 2:nrow(Xnmp)){
				(year <- Xnmp$yr[i])
				(year_m1 <- ifelse(i==1,Xnmp$yr[1],Xnmp$yr[i]-1))
			(sub <- subset(Xnmp, yr==year))
			(sub_m1 <- subset(Xnmp, yr==year_m1))
		(Xnmp$HypTr[i] <- ifelse(sub$trade_withNonMbrPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	Xnmp[,c("yr","trade_withNonMbrPpts","HypTr")]
	dm2 <- Xnmp[,c("yr","trade_withNonMbrPpts","HypTr")];names(dm2)[3]<-"hyptr_xdata_NMP"; head(dm2)

	#Export data, trade with non-ppts
	(Xnp <- merge(NP,subNP, by.x="yr", by.y="yrM1", all.x=T))
			Xnp$tradeincrease <- Xnp[,2] - Xnp[,3]; Xnp
			Xnp$tradeincrease[1] <- 0; Xnp

			Xnp$adj <- rep(0,nrow(Xnp))					#hyp 		/ obs
		Xnp$adj[Xnp$yr<=1988] 				<- 1.00 / 1.00	#(NP-NP) 	 / (NP-NP)
		Xnp$adj[Xnp$yr>=1989 & Xnp$yr<=2000] 	<- 1.22 / 1.00	#(**M**-NP) / (NP-NP)
		Xnp$adj[Xnp$yr>=2001] 				<- 1.22 / 1.22	#(M-NP)	 / ( M-NP)
			Xnp$adjustment <- rep(0,nrow(Xnp))
		Xnp$adjustment[Xnp$yr<=1988] 				<- "=1.00/1.00"
		Xnp$adjustment[Xnp$yr>=1989 & Xnp$yr<=2000] 	<- "=1.22/1.00"
		Xnp$adjustment[Xnp$yr>=2001] 				<- "=1.22/1.22"

		Xnp$HypTr <- rep(0,nrow(Xnp))
			Xnp$HypTr[1] <- Xnp$trade_withNonPpts[1]
		for(i in 2:nrow(Xnp)){
				(year <- Xnp$yr[i])
				(year_m1 <- ifelse(i==1,Xnp$yr[1],Xnp$yr[i]-1))
			(sub <- subset(Xnp, yr==year))
			(sub_m1 <- subset(Xnp, yr==year_m1))
		(Xnp$HypTr[i] <- ifelse(sub$trade_withNonPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	Xnp[,c("yr","trade_withNonPpts","HypTr")]
	dm3 <- Xnp[,c("yr","trade_withNonPpts","HypTr")];names(dm3)[3]<-"hyptr_xdata_NP"; head(dm3)

	#Country imports
	ctyM <- subset(grt, tomz_1 == 621)
		head(ctyM)
	(ctyMyrs<-Myrs <- as.data.frame(sort(unique(ctyM$year)))) 	#1950-2003
	names(Myrs) <- "yr"
	Myrs$trade_withMbrs <-rep(0, length(ctyMyrs))
	Myrs$trade_withNonMbrPpts<-rep(0, length(ctyMyrs))
	Myrs$trade_withNonPpts<-rep(0, length(ctyMyrs))
		head(Myrs)
	for(i in 1:nrow(Myrs)){
		sub <- subset(ctyM, year==Myrs$yr[i])
		Myrs$trade_withMbrs[i] <- ifelse(sum(sub$mbr_ind_2==1)==0,0,sum(subset(sub,mbr_ind_2==1)$imports_nonlog))
		Myrs$trade_withNonMbrPpts[i] <- ifelse(sum(sub$nonmbrppt_ind_2==1)==0,0,sum(subset(sub,nonmbrppt_ind_2==1)$imports_nonlog))
		Myrs$trade_withNonPpts[i] <- ifelse(sum(sub$nonppt_ind_2==1)==0,0,sum(subset(sub,nonppt_ind_2==1)$imports_nonlog))
		}
	head(Myrs)
		M 	<- Myrs[,c("yr","trade_withMbrs")]
		NMP 	<- Myrs[,c("yr","trade_withNonMbrPpts")]
		NP	<- Myrs[,c("yr","trade_withNonPpts")]

		M1 <- Myrs
	M1$yrM1 <- M1$yr+1		
			head(M1)
		subM <- M1[,c("yrM1","trade_withMbrs")]; names(subM)[2]<-"twM_m1"; head(subM)
		subNMP <- M1[,c("yrM1","trade_withNonMbrPpts")]; names(subNMP)[2]<-"twNMP_m1"; head(subNMP)
		subNP <- M1[,c("yrM1","trade_withNonPpts")]; names(subNP)[2]<-"twNP_m1"; head(subNP)

	#Import data, trade with formal members
	(Mm <- merge(M, subM, by.x="yr", by.y="yrM1", all.x=T))
			Mm$tradeincrease <- Mm[,2] - Mm[,3]; Mm
			Mm$tradeincrease[1] <- 0; Mm

			Mm$adj <- rep(0,nrow(Mm))				#hyp / obs
		Mm$adj[Mm$yr<=1988] 			<- 1.22 / 1.22	#(NP-M) 	/ (NP-M)
		Mm$adj[Mm$yr>=1989 & Mm$yr<=2000] 	<- 1.41 / 1.22	#(**M**-M) 	/ (NP-M)
		Mm$adj[Mm$yr>=2001] 			<- 1.41 / 1.41	#(M-M)	/ (M-M)
			Mm$adjustment <- rep(0,nrow(Mm))
		Mm$adjustment[Mm$yr<=1988] 			<- "=1.22/1.22"
		Mm$adjustment[Mm$yr>=1989 & Mm$yr<=2000] 	<- "=1.41/1.22"
		Mm$adjustment[Mm$yr>=2001] 			<- "=1.41/1.41"

		Mm$HypTr <- rep(0,nrow(Mm))
			Mm$HypTr[1] <- Mm$trade_withMbrs[1]
		for(i in 2:nrow(Mm)){
				(year <- Mm$yr[i])
				(year_m1 <- ifelse(i==1,Mm$yr[1],Mm$yr[i]-1))
			(sub <- subset(Mm, yr==year))
			(sub_m1 <- subset(Mm, yr==year_m1))
		(Mm$HypTr[i] <- ifelse(sub$trade_withMbrs==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	dm4 <- Mm[,c("yr","trade_withMbrs","HypTr")];names(dm4)[3]<-"hyptr_mdata_M"; head(dm4)

	#Import data, trade with non-mbr ppts
	(Mnmp <- merge(NMP,subNMP, by.x="yr", by.y="yrM1", all.x=T))
			Mnmp$tradeincrease <- Mnmp[,2] - Mnmp[,3]; Mnmp
			Mnmp$tradeincrease[1] <- 0; Mnmp

			Mnmp$adj <- rep(0,nrow(Mnmp))					#hyp 		/ obs
		Mnmp$adj[Mnmp$yr<=1988] 			<- 1.19 / 1.19	#(NP-NMP) 	 / (NP-NMP)
		Mnmp$adj[Mnmp$yr>=1989 & Mnmp$yr<=2000] 	<- 1.46 / 1.19	#(**M**-NMP) / (NP-NMP)
		Mnmp$adj[Mnmp$yr>=2001] 			<- 1.46 / 1.46	#(M-NMP)	 / ( M-NMP)
			Mnmp$adjustment <- rep(0,nrow(Mnmp))
		Mnmp$adjustment[Mnmp$yr<=1988] 			<- "=1.19/1.19"
		Mnmp$adjustment[Mnmp$yr>=1989 & Mnmp$yr<=2000] 	<- "=1.46/1.19"
		Mnmp$adjustment[Mnmp$yr>=2001] 			<- "=1.46/1.46"

		Mnmp$HypTr <- rep(0,nrow(Mnmp))
			Mnmp$HypTr[1] <- Mnmp$trade_withNonMbrPpts[1]
		for(i in 2:nrow(Mnmp)){
				(year <- Mnmp$yr[i])
				(year_m1 <- ifelse(i==1,Mnmp$yr[1],Mnmp$yr[i]-1))
			(sub <- subset(Mnmp, yr==year))
			(sub_m1 <- subset(Mnmp, yr==year_m1))
		(Mnmp$HypTr[i] <- ifelse(sub$trade_withNonMbrPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	Mnmp[,c("yr","trade_withNonMbrPpts","HypTr")]
	dm5 <- Mnmp[,c("yr","trade_withNonMbrPpts","HypTr")];names(dm5)[3]<-"hyptr_mdata_NMP"; head(dm5)

	#Export data, trade with non-ppts
	(Mnp <- merge(NP,subNP, by.x="yr", by.y="yrM1", all.x=T))
			Mnp$tradeincrease <- Mnp[,2] - Mnp[,3]; Mnp
			Mnp$tradeincrease[1] <- 0; Mnp

			Mnp$adj <- rep(0,nrow(Mnp))					#hyp 		/ obs
		Mnp$adj[Mnp$yr<=1988] 				<- 1.00 / 1.00	#(NP-NP) 	 / (NP-NP)
		Mnp$adj[Mnp$yr>=1989 & Mnp$yr<=2000] 	<- 1.22 / 1.00	#(**M**-NP) / (NP-NP)
		Mnp$adj[Mnp$yr>=2001] 				<- 1.22 / 1.22	#(M-NP)	 / ( M-NP)
			Mnp$adjustment <- rep(0,nrow(Mnp))
		Mnp$adjustment[Mnp$yr<=1988] 				<- "=1.00/1.00"
		Mnp$adjustment[Mnp$yr>=1989 & Mnp$yr<=2000] 	<- "=1.22/1.00"
		Mnp$adjustment[Mnp$yr>=2001] 				<- "=1.22/1.22"

		Mnp$HypTr <- rep(0,nrow(Mnp))
			Mnp$HypTr[1] <- Mnp$trade_withNonPpts[1]
		for(i in 2:nrow(Mnp)){
				(year <- Mnp$yr[i])
				(year_m1 <- ifelse(i==1,Mnp$yr[1],Mnp$yr[i]-1))
			(sub <- subset(Mnp, yr==year))
			(sub_m1 <- subset(Mnp, yr==year_m1))
		(Mnp$HypTr[i] <- ifelse(sub$trade_withNonPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	Mnp[,c("yr","trade_withNonPpts","HypTr")]
	dm6 <- Mnp[,c("yr","trade_withNonPpts","HypTr")];names(dm6)[3]<-"hyptr_mdata_NP"; head(dm6)
	dm6

	(yrs <- sort(unique(grt$year))) #1946-2004
	h <- as.data.frame(yrs)
	head(h)
	head(dm1)
	(h1 <- merge(h, dm1, by.x="yrs",by.y="yr",all.x=T))
	(h2 <- merge(h1, dm2, by.x="yrs",by.y="yr",all.x=T))
	(h3 <- merge(h2, dm3, by.x="yrs",by.y="yr",all.x=T))
	(h4 <- merge(h3, dm4, by.x="yrs",by.y="yr",all.x=T))
	(h5 <- merge(h4, dm5, by.x="yrs",by.y="yr",all.x=T))
	(h6 <- merge(h5, dm6, by.x="yrs",by.y="yr",all.x=T))

	head(h6)
	h6 <- as.data.frame(h6)
	h6$obs <- rep(NA,nrow(h6))
	for(i in 1:nrow(h6)){
	h6$obs[i] <- sum(h6$trade_withMbrs.x[i],
				h6$trade_withMbrs.y[i],
			h6$trade_withNonMbrPpts.x[i], 
				h6$trade_withNonMbrPpts.y[i], 
			h6$trade_withNonPpts.x[i], 
				h6$trade_withNonPpts.y[i], na.rm=T)
			}
		h6[,c("yrs","obs")]
	h6$hyp <- rep(NA,nrow(h6))
	for(i in 1:nrow(h6)){
	h6$hyp[i] <- sum(h6$hyptr_xdata_M[i],
				h6$hyptr_mdata_M[i],
			h6$hyptr_xdata_NMP[i], 
				h6$hyptr_mdata_NMP[i], 
			h6$hyptr_xdata_NP[i], 
				h6$hyptr_mdata_NP[i], na.rm=T)
			}
		h6[,c("yrs","hyp")]
	h6$obs[h6$yrs==2004 | h6$yrs<=1949] <- h6$hyp[h6$yrs==2004 | h6$yrs<=1949] <- NA
	h6[,c("yrs","obs","hyp")]
#pdf("JOP_Appendix_FigureA1_Left_China.pdf")
	plot(h6$yrs,h6$hyp,type="n",axes=F,xlab="",ylab="",ylim=c(0,500000000000),xlim=c(1945,2015))
		axis(1); axis(2,at=c(0, 100000000000, 200000000000,
						300000000000, 400000000000, 500000000000),
				c("0","$100","$200","$300","$400","$500"),las=2)
		points(h6$yrs,h6$obs,font=2,pch=16)
		points(h6$yrs,h6$hyp,font=2,pch=3,lwd=2)
		mtext("year",font=2,line=2.5,side=1)
		mtext("China",font=2,cex=1.75,line=2.5,side=3)
		mtext("trade ($billions)",font=2,line=1,side=3,cex=1.25)
			segments(1988.5, 0, 	1988.5, 500000000000,lty=3)
			segments(2001.5, 0, 	2001.5, 500000000000,lty=3)
		text(2004, 180000000000, "observed",font=2,adj=c(0,0),cex=.85)
		text(2004, 470000000000, "hypothetical",font=2,adj=c(0,0),cex=.85)
#dev.off()

	#China
	h6[,c("yrs","obs","hyp")]
h6[c(43,56),c("yrs","obs","hyp")]

	#Observed values @ years={1988, 2001}
	h6[43,]$obs	#o1 = observed values in 1988
o1 <-  32155365808	# 32 155 365 808	$32.2 bill
	h6[56,]$obs	#o1 = observed values in 2001
o2 <- 115241846849 	#115 241 846 849 	$115.2 bill

	#Hypothetical increase, 1988 -- 2001:
	h6[43,]$hyp		#h1 = hypothetical value in 1988
h1 <-  32155365808 	# 32 155 365 808	$ 32.2 bill
	h6[56,]$hyp		#h1 = hypothetical value in 2001
h2 <- 389631585951 	#389 631 585 951	$389.6 bill

	(o2-o1)/o1				#2.58	#grew 258% over this period
	(o2-o1)/o1/(2001-1988+1)	#simple average growth of 18.5% per year. . . 

	(h2-h1)/h1				#11.117 #potentially grow 1281%
	(h2-h1)/h1/(2001-1988+1)	#simple average growth of 79.4% per year. . . 

	#Comparison, observed versus hypothetical trade, 2001. . . (LR effect)
(h2-o2)/o2	#==238.1% boost in total trade in 2001. . . 


###AGGREGATE TRADE
#MEXICO, 1979 versus 1986. . .

	(yrs <- sort(unique(grt$year))) #1946-2004

	ctyX <- subset(grt, tomz_2 == 411)
	head(ctyX)
		###1946--85: NP: "out"
		###1986--04: M : "art33", "wto"

		unique(ctyX[,c("year","gatt_2")])
	(Xyrs <- as.data.frame(sort(unique(ctyX$year)))) 	#1946-2004
	names(Xyrs) <- "yr"
	Xyrs$trade_withMbrs <-rep(0, nrow(Xyrs))
	Xyrs$trade_withNonMbrPpts<-rep(0, nrow(Xyrs))
	Xyrs$trade_withNonPpts<-rep(0, nrow(Xyrs))
		head(Xyrs)
	for(i in 1:nrow(Xyrs)){
		sub <- subset(ctyX, year==Xyrs$yr[i])
		Xyrs$trade_withMbrs[i] <- ifelse(sum(sub$mbr_ind_1==1)==0,0,sum(subset(sub,mbr_ind_1==1)$imports_nonlog))
		Xyrs$trade_withNonMbrPpts[i] <- ifelse(sum(sub$nonmbrppt_ind_1==1)==0,0,sum(subset(sub,nonmbrppt_ind_1==1)$imports_nonlog))
		Xyrs$trade_withNonPpts[i] <- ifelse(sum(sub$nonppt_ind_1==1)==0,0,sum(subset(sub,nonppt_ind_1==1)$imports_nonlog))
		}
	head(Xyrs)
		M 	<- Xyrs[,c("yr","trade_withMbrs")]
		NMP 	<- Xyrs[,c("yr","trade_withNonMbrPpts")]
		NP	<- Xyrs[,c("yr","trade_withNonPpts")]

		M1 <- Xyrs
	M1$yrM1 <- M1$yr+1		
			head(M1)
		subM <- M1[,c("yrM1","trade_withMbrs")]; names(subM)[2]<-"twM_m1"; head(subM)
		subNMP <- M1[,c("yrM1","trade_withNonMbrPpts")]; names(subNMP)[2]<-"twNMP_m1"; head(subNMP)
		subNP <- M1[,c("yrM1","trade_withNonPpts")]; names(subNP)[2]<-"twNP_m1"; head(subNP)

	#Export data, trade with formal members
	(Xm <- merge(M, subM, by.x="yr", by.y="yrM1", all.x=T))
			Xm$tradeincrease <- Xm[,2] - Xm[,3]; Xm
			Xm$tradeincrease[1] <- 0; Xm

			Xm$adj <- rep(0,nrow(Xm))				#hyp / obs
		Xm$adj[Xm$yr<=1978] 			<- 1.0 / 1.0	#(NP-NP) / (NP-NP)
		Xm$adj[Xm$yr>=1979 & Xm$yr<=1985] 	<- 1.41 / 1.22	#(**M**-M) 	/ (NP-M)
		Xm$adj[Xm$yr>=1986] 			<- 1.41 / 1.41	#(M-M)	/ (M-M)
			Xm$adjustment <- rep(0,nrow(Xm))
		Xm$adjustment[Xm$yr<=1978] 			<- "=1.00/1.00"
		Xm$adjustment[Xm$yr>=1979 & Xm$yr<=1985] 	<- "=1.41/1.22"
		Xm$adjustment[Xm$yr>=1986] 			<- "=1.41/1.41"

			Xm
		(Xm <- subset(Xm,yr>=1948)); sort(unique(Xm$yr))
		Xm$HypTr <- rep(0,nrow(Xm))
			Xm$HypTr[1] <- Xm$trade_withMbrs[1]
				Xm
		for(i in 2:nrow(Xm)){
				(year <- Xm$yr[i])
				(year_m1 <- ifelse(i==1,Xm$yr[1],Xm$yr[i]-1))
			(sub <- subset(Xm, yr==year))
			(sub_m1 <- subset(Xm, yr==year_m1))
		(Xm$HypTr[i] <- ifelse(sub$trade_withMbrs==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	dm1 <- Xm[,c("yr","trade_withMbrs","HypTr")];names(dm1)[3]<-"hyptr_xdata_M"; head(dm1)

	#Export data, trade with non-mbr ppts
	(Xnmp <- merge(NMP,subNMP, by.x="yr", by.y="yrM1", all.x=T))
			Xnmp$tradeincrease <- Xnmp[,2] - Xnmp[,3]; Xnmp
			Xnmp$tradeincrease[1] <- 0; Xnmp

			Xnmp$adj <- rep(0,nrow(Xnmp))				#hyp / obs
		Xnmp$adj[Xnmp$yr<=1978] 			<- 1.19 / 1.19	#(NP-NMP) / (NP-NMP)
		Xnmp$adj[Xnmp$yr>=1979 & Xnmp$yr<=1985] 	<- 1.46 / 1.19	#(**M**-NMP) 	/ (NP-NMP)
		Xnmp$adj[Xnmp$yr>=1986] 			<- 1.46 / 1.46	#(M-NMP)	/ (M-NMP)
			Xnmp$adjustment <- rep(0,nrow(Xnmp))
		Xnmp$adjustment[Xnmp$yr<=1978] 			<- "=1.19/1.19"
		Xnmp$adjustment[Xnmp$yr>=1979 & Xnmp$yr<=1985] 	<- "=1.46/1.19"
		Xnmp$adjustment[Xnmp$yr>=1986] 			<- "=1.46/1.46"

		Xnmp$HypTr <- rep(0,nrow(Xnmp))
			Xnmp
		(Xnmp <- subset(Xnmp,yr>=1948 & yr<=1996)); sort(unique(Xnmp$yr))
			Xnmp$HypTr[1] <- Xnmp$trade_withNonMbrPpts[1]
		for(i in 2:nrow(Xnmp)){
				(year <- Xnmp$yr[i])
				(year_m1 <- ifelse(i==1,Xnmp$yr[1],Xnmp$yr[i]-1))
			(sub <- subset(Xnmp, yr==year))
			(sub_m1 <- subset(Xnmp, yr==year_m1))
		(Xnmp$HypTr[i] <- ifelse(sub$trade_withNonMbrPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
Xnmp
	Xnmp[,c("yr","trade_withNonMbrPpts","HypTr")]
	dm2 <- Xnmp[,c("yr","trade_withNonMbrPpts","HypTr")];names(dm2)[3]<-"hyptr_xdata_NMP"; head(dm2)

	#Export data, trade with non-ppts
	(Xnp <- merge(NP,subNP, by.x="yr", by.y="yrM1", all.x=T))
			Xnp$tradeincrease <- Xnp[,2] - Xnp[,3]; Xnp
			Xnp$tradeincrease[1] <- 0; Xnp

			Xnp$adj <- rep(0,nrow(Xnp))					#hyp 		/ obs
		Xnp$adj[Xnp$yr<=1978] 				<- 1.00 / 1.00	#(NP-NP) 	 / (NP-NP)
		Xnp$adj[Xnp$yr>=1979 & Xnp$yr<=1985] 	<- 1.22 / 1.00	#(**M**-NP) / (NP-NP)
		Xnp$adj[Xnp$yr>=1986] 				<- 1.22 / 1.22	#(M-NP)	 / ( M-NP)
			Xnp$adjustment <- rep(0,nrow(Xnp))
		Xnp$adjustment[Xnp$yr<=1978] 				<- "=1.00/1.00"
		Xnp$adjustment[Xnp$yr>=1979 & Xnp$yr<=1985] 	<- "=1.22/1.00"
		Xnp$adjustment[Xnp$yr>=1986] 				<- "=1.22/1.22"

		Xnp$HypTr <- rep(0,nrow(Xnp))
			Xnp$HypTr[1] <- Xnp$trade_withNonPpts[1]
Xnp
		for(i in 2:nrow(Xnp)){
				(year <- Xnp$yr[i])
				(year_m1 <- ifelse(i==1,Xnp$yr[1],Xnp$yr[i]-1))
			(sub <- subset(Xnp, yr==year))
			(sub_m1 <- subset(Xnp, yr==year_m1))
		(Xnp$HypTr[i] <- ifelse(sub$trade_withNonPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	Xnp[,c("yr","trade_withNonPpts","HypTr")]
	dm3 <- Xnp[,c("yr","trade_withNonPpts","HypTr")];names(dm3)[3]<-"hyptr_xdata_NP"; head(dm3)

	ctyM <- subset(grt, tomz_1 == 411)
		head(ctyM)
	(Myrs <- as.data.frame(sort(unique(ctyM$year)))) 	#1950-2003
	names(Myrs) <- "yr"
	Myrs$trade_withMbrs <-rep(0,nrow(Myrs))
	Myrs$trade_withNonMbrPpts<-rep(0,nrow(Myrs))
	Myrs$trade_withNonPpts<-rep(0,nrow(Myrs))
		head(Myrs)
	for(i in 1:nrow(Myrs)){
		sub <- subset(ctyM, year==Myrs$yr[i])
		Myrs$trade_withMbrs[i] <- ifelse(sum(sub$mbr_ind_2==1)==0,0,sum(subset(sub,mbr_ind_2==1)$imports_nonlog))
		Myrs$trade_withNonMbrPpts[i] <- ifelse(sum(sub$nonmbrppt_ind_2==1)==0,0,sum(subset(sub,nonmbrppt_ind_2==1)$imports_nonlog))
		Myrs$trade_withNonPpts[i] <- ifelse(sum(sub$nonppt_ind_2==1)==0,0,sum(subset(sub,nonppt_ind_2==1)$imports_nonlog))
		}
	head(Myrs)
		M 	<- Myrs[,c("yr","trade_withMbrs")]
		NMP 	<- Myrs[,c("yr","trade_withNonMbrPpts")]
		NP	<- Myrs[,c("yr","trade_withNonPpts")]

		M1 <- Myrs
	M1$yrM1 <- M1$yr+1		
			head(M1)
		subM <- M1[,c("yrM1","trade_withMbrs")]; names(subM)[2]<-"twM_m1"; head(subM)
		subNMP <- M1[,c("yrM1","trade_withNonMbrPpts")]; names(subNMP)[2]<-"twNMP_m1"; head(subNMP)
		subNP <- M1[,c("yrM1","trade_withNonPpts")]; names(subNP)[2]<-"twNP_m1"; head(subNP)

	#Import data, trade with formal members
	(Mm <- merge(M, subM, by.x="yr", by.y="yrM1", all.x=T))
			Mm$tradeincrease <- Mm[,2] - Mm[,3]; Mm
			Mm$tradeincrease[1] <- 0; Mm

			Mm$adj <- rep(0,nrow(Mm))				#hyp / obs
		Mm$adj[Mm$yr<=1978] 			<- 1.22 / 1.22	#(NP-M) 	/ (NP-M)
		Mm$adj[Mm$yr>=1979 & Mm$yr<=1985] 	<- 1.41 / 1.22	#(**M**-M) 	/ (NP-M)
		Mm$adj[Mm$yr>=1986] 			<- 1.41 / 1.41	#(M-M)	/ (M-M)
			Mm$adjustment <- rep(0,nrow(Mm))
		Mm$adjustment[Mm$yr<=1978] 			<- "=1.22/1.22"
		Mm$adjustment[Mm$yr>=1979 & Mm$yr<=1985] 	<- "=1.41/1.22"
		Mm$adjustment[Mm$yr>=1986] 			<- "=1.41/1.41"

		Mm$HypTr <- rep(0,nrow(Mm))
			Mm
			(Mm <-subset(Mm, yr >=1948))
			Mm$HypTr[1] <- Mm$trade_withMbrs[1]
		for(i in 2:nrow(Mm)){
				(year <- Mm$yr[i])
				(year_m1 <- ifelse(i==1,Mm$yr[1],Mm$yr[i]-1))
			(sub <- subset(Mm, yr==year))
			(sub_m1 <- subset(Mm, yr==year_m1))
		(Mm$HypTr[i] <- ifelse(sub$trade_withMbrs==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	dm4 <- Mm[,c("yr","trade_withMbrs","HypTr")];names(dm4)[3]<-"hyptr_mdata_M"; head(dm4)

	#Import data, trade with non-mbr ppts
	(Mnmp <- merge(NMP,subNMP, by.x="yr", by.y="yrM1", all.x=T))
			Mnmp$tradeincrease <- Mnmp[,2] - Mnmp[,3]; Mnmp
			Mnmp$tradeincrease[1] <- 0; Mnmp

			Mnmp$adj <- rep(0,nrow(Mnmp))					#hyp 		/ obs
		Mnmp$adj[Mnmp$yr<=1978] 			<- 1.19 / 1.19	#(NP-NMP) 	 / (NP-NMP)
		Mnmp$adj[Mnmp$yr>=1979 & Mnmp$yr<=1985] 	<- 1.46 / 1.19	#(**M**-NMP) / (NP-NMP)
		Mnmp$adj[Mnmp$yr>=1986] 			<- 1.46 / 1.46	#(M-NMP)	 / ( M-NMP)
			Mnmp$adjustment <- rep(0,nrow(Mnmp))
		Mnmp$adjustment[Mnmp$yr<=1978] 			<- "=1.19/1.19"
		Mnmp$adjustment[Mnmp$yr>=1979 & Mnmp$yr<=1985] 	<- "=1.46/1.19"
		Mnmp$adjustment[Mnmp$yr>=1986] 			<- "=1.46/1.46"

		Mnmp$HypTr <- rep(0,nrow(Mnmp))
	Mnmp
			Mnmp$HypTr[1] <- Mnmp$trade_withNonMbrPpts[1]
		for(i in 2:nrow(Mnmp)){
				(year <- Mnmp$yr[i])
				(year_m1 <- ifelse(i==1,Mnmp$yr[1],Mnmp$yr[i]-1))
			(sub <- subset(Mnmp, yr==year))
			(sub_m1 <- subset(Mnmp, yr==year_m1))
		(Mnmp$HypTr[i] <- ifelse(sub$trade_withNonMbrPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	Mnmp[,c("yr","trade_withNonMbrPpts","HypTr")]
	dm5 <- Mnmp[,c("yr","trade_withNonMbrPpts","HypTr")];names(dm5)[3]<-"hyptr_mdata_NMP"; head(dm5)

	#Export data, trade with non-ppts
	(Mnp <- merge(NP,subNP, by.x="yr", by.y="yrM1", all.x=T))
			Mnp$tradeincrease <- Mnp[,2] - Mnp[,3]; Mnp
			Mnp$tradeincrease[1] <- 0; Mnp

			Mnp$adj <- rep(0,nrow(Mnp))					#hyp 		/ obs
		Mnp$adj[Mnp$yr<=1978] 				<- 1.00 / 1.00	#(NP-NP) 	 / (NP-NP)
		Mnp$adj[Mnp$yr>=1979 & Mnp$yr<=1985] 	<- 1.22 / 1.00	#(**M**-NP) / (NP-NP)
		Mnp$adj[Mnp$yr>=1986] 				<- 1.22 / 1.22	#(M-NP)	 / ( M-NP)
			Mnp$adjustment <- rep(0,nrow(Mnp))
		Mnp$adjustment[Mnp$yr<=1978] 				<- "=1.00/1.00"
		Mnp$adjustment[Mnp$yr>=1979 & Mnp$yr<=1985] 	<- "=1.22/1.00"
		Mnp$adjustment[Mnp$yr>=1986] 				<- "=1.22/1.22"

		Mnp$HypTr <- rep(0,nrow(Mnp))
	Mnp
			Mnp$HypTr[1] <- Mnp$trade_withNonPpts[1]
		for(i in 2:nrow(Mnp)){
				(year <- Mnp$yr[i])
				(year_m1 <- ifelse(i==1,Mnp$yr[1],Mnp$yr[i]-1))
			(sub <- subset(Mnp, yr==year))
			(sub_m1 <- subset(Mnp, yr==year_m1))
		(Mnp$HypTr[i] <- ifelse(sub$trade_withNonPpts==0,0,(sub_m1$HypTr[1]+ sub$tradeincrease[1])*sub$adj[1]))
		}
	Mnp[,c("yr","trade_withNonPpts","HypTr")]
	dm6 <- Mnp[,c("yr","trade_withNonPpts","HypTr")];names(dm6)[3]<-"hyptr_mdata_NP"; head(dm6)
dm6

	(yrs <- sort(unique(grt$year))) #1946-2004
	h <- as.data.frame(yrs)
	head(h)
	head(dm1)
	(h1 <- merge(h, dm1, by.x="yrs",by.y="yr",all.x=T))
	(h2 <- merge(h1, dm2, by.x="yrs",by.y="yr",all.x=T))
	(h3 <- merge(h2, dm3, by.x="yrs",by.y="yr",all.x=T))
	(h4 <- merge(h3, dm4, by.x="yrs",by.y="yr",all.x=T))
	(h5 <- merge(h4, dm5, by.x="yrs",by.y="yr",all.x=T))
	(h6 <- merge(h5, dm6, by.x="yrs",by.y="yr",all.x=T))

	head(h6)
	h6 <- as.data.frame(h6)
	h6$obs <- rep(NA,nrow(h6))
	for(i in 1:nrow(h6)){
	h6$obs[i] <- sum(h6$trade_withMbrs.x[i],
				h6$trade_withMbrs.y[i],
			h6$trade_withNonMbrPpts.x[i], 
				h6$trade_withNonMbrPpts.y[i], 
			h6$trade_withNonPpts.x[i], 
				h6$trade_withNonPpts.y[i], na.rm=T)
			}
		h6[,c("yrs","obs")]
	h6$hyp <- rep(NA,nrow(h6))
	for(i in 1:nrow(h6)){
	h6$hyp[i] <- sum(h6$hyptr_xdata_M[i],
				h6$hyptr_mdata_M[i],
			h6$hyptr_xdata_NMP[i], 
				h6$hyptr_mdata_NMP[i], 
			h6$hyptr_xdata_NP[i], 
				h6$hyptr_mdata_NP[i], na.rm=T)
			}
		h6[,c("yrs","hyp")]
	#h6$obs[h6$yrs==2004 | h6$yrs<=1949] <- h6$hyp[h6$yrs==2004 | h6$yrs<=1949] <- NA
	h6[,c("yrs","obs","hyp")]

#pdf("JOP_Appendix_FigureA1_Right_Mexico.pdf")
	plot(h6$yrs,h6$hyp,type="n",axes=F,xlab="",ylab="",ylim=c(0,100000000000),xlim=c(1945,2015))
		axis(1); axis(2,at=c(0, 10000000000, 20000000000,
						30000000000, 40000000000, 50000000000,
						60000000000, 70000000000,
						80000000000, 90000000000, 100000000000),
				c("0","$10","$20","$30","$40","$50","$60","$70","$80","$90","$100"),las=2)
		points(h6$yrs,h6$obs,font=2,pch=16)
		points(h6$yrs,h6$hyp,font=2,pch=3,lwd=2)
		mtext("year",font=2,line=2.5,side=1)
		mtext("Mexico",font=2,cex=1.75,line=2.5,side=3)
		mtext("trade ($billions)",font=2,line=1,side=3,cex=1.25)
			segments(1978.5, 0, 	1978.5, 500000000000,lty=3)
			segments(1985.5, 0, 	1985.5, 500000000000,lty=3)
		text(2002, 55000000000, "observed",font=2,adj=c(0,0),cex=.85)
		text(2002, 95000000000, "hypothetical",font=2,adj=c(0,0),cex=.85)
#dev.off()

	#Mexico
	h6[,c("yrs","obs","hyp")]

h6[c(33,41),c("yrs","obs","hyp")]
	#Observed increase, 1978 -- 1986:
o1 <-  8134019381	# 8 134 019 381 == $8.1 bill
o2 <- 11075115269		# 11075 115 269 	== $11.1 bill

	#Hypothetical increase, 1978 -- 1986:
h1 <-  8134019381	# 8134019381 == $ 8.1 bill
h2 <- 35027380264	#35 027380 264== $35.0 bill

	(o2-o1)/o1				#.36	#grew 36%
	(o2-o1)/o1/(1986-1978+1)	#simple average growth of 4% per year. . . 

	(h2-h1)/h1				#3.31 #potentially grow 331%
	(h2-h1)/h1/(1986-1978+1)	#simple average growth of 37% per year. . . 

	#Comparison, observed versus hypothetical trade @1986. . . (LR effect)
o2 <- 11075115269		# 11075 115 269 	== $11.1 bill
h2 <- 35027380264	#35 027380 264== $35.0 bill

(h2-o2)/o2	#==2.1627% boost in total trade in 2004. . . 

	#Comparison, observed versus hypothetical trade, 2004. . . (LR effect)
		#Observed increase, 1978 -- 2004:
	o2 <- 60697006232	# 60 697 006 232	#60.697 bill
		#Hypothetical increase, 1978 -- 2004:
	h2 <- 84564445460 #84 564 445 460	#84.56 Bill 
(h2-o2)/o2	#==39% boost in total trade in 2004. . . 


#######
# Figure A.2
#	World Membership Percent, 1948-2015

	years <- 1948:2015
	percMbr <- rep(NA,length(years))
	for(i in 1:length(years)){
	subsub <- subset(d10, year==years[i] & cty_ind==1)
	(MbrCt <- nrow(subset(subsub,JoinYear_PrevJoined==1 | JoinYear_Ind==1)))
	(TotCtys    <- nrow(subsub))
	percMbr[i] <- MbrCt / (TotCtys)
	}

# pdf("JOP_Appendix_FigureA2_PercWldMbrs.pdf")
		plot(years,percMbr, type="l",lwd=2,ylim=c(0,1), xlim=c(1944,2016),axes=F,
				xlab="",ylab="")
		axis(1,seq(1945,2015,10))
		axis(2, at=seq(0, 1, .2), c("0%","20%","40%","60%","80%","100%"), las=2)
		abline(v=1994.5, lty=3)
		text(1994, .8,"GATT", adj=c(1,1),font=2)
		text(1995, .8,"WTO",adj=c(0,1),font=2)
		mtext("Year", side=1, line=2.5)
		#mtext("Percent Member Countries", side=2, line=2.5)
		mtext("Percent World Members", side=3, line=1,font=2,cex=1.5)
#	dev.off()

#######
# Figure A.3
#	Geopolitical Alignment and Accession Rigor
#	Using data from Allee and Scalera (2012)

alleescalera <- read.csv("AlleeScalera_IO_2012.csv")

	(rigacc <- subset(alleescalera, rigorous==1 & is.na(year)==T)[,c("country",
					"cowcode","rigaccyears",
					"rignumwpmems_max",
					"rigdiffappacc1",
					"rigcommit_pg")])

	rigacc$cowcode[rigacc$country=="Tonga"] <- 955
	rigacc$cowcode[rigacc$country=="Taiwan,China"] <- 713
	rigacc$AlleeScalera_Ind <- rep(1, nrow(rigacc))
		names(rigacc)
	(rigacc_merge <- rigacc[,c("cowcode",
					"rignumwpmems_max",
					"rigdiffappacc1",
					"rigcommit_pg","AlleeScalera_Ind")])
	(names(rigacc_merge) <- c("ccode",
					"AS_rignumwpmems_max",
					"AS_rigdiffappacc1",
					"AS_rigcommit_pg","AlleeScalera_Ind"))

	d10merge2 <- merge(d10merge,rigacc_merge,by="ccode",all.x=T)
		names(d10merge2)

	h <- subset(d10merge2, (ApplyYrI_Ind_a==1 & year>=1995 & AlleeScalera_Ind==1) |
		(ApplyYrI_Ind_a==1 &AccessionType=="GATT, Art33" & AlleeScalera_Ind==1)|
		(ApplyYrI_Ind_a==1 &AccessionType=="WTO, Art12" & AlleeScalera_Ind==1))

	nrow(h)	#48 countries have Allee and Scalera data

#pdf("JOP_Appendix_FigureA3_Left_AppliedTariffs.pdf")
	plot(h$s3un100, h$AS_rigdiffappacc1, type="n",xlim=c(-100,100),ylim=c(-10,35),
		xlab="",ylab="",axes=F); axis(1); axis(2,las=2)
		#text(jitter(h$s3un100), jitter(h$AS_rigdiffappacc1), h$stateabb,cex=.75,font=2)
		text(jitter(h$s3un100), jitter(h$AS_rigdiffappacc1), h$statenme,cex=.75,font=2)
	summary(fit <- lm(AS_rigdiffappacc1 ~ s3un100, data=h))	
		abline(summary(fit))
		text(90, -10,"(-) 10%", font=2,cex=1.5)
	axis(1);axis(2,las=2)
	mtext("UN Voting Similarity",side=1,line=2.5,font=2,cex=1.25)	
	mtext("Applied Tariff Decrease",side=3,line=2.5,font=2,cex=1.5)	
	mtext("(Allee Scalera (2011) data)",side=3,line=1,font=2,cex=1)	
		nobs <-as.data.frame(cbind(h$s3un100,h$AS_rigdiffappacc1))
		nrow(subset(nobs, is.na(nobs[,1])==F & is.na(nobs[,2])==F))
		#45 observations
	mtext("(n=45)",side=3,line=0,font=2,cex=1)
#dev.off()

#pdf("JOP_Appendix_FigureA3_Center_WpMbrCt.pdf")
	plot(h$s3un100, h$AS_rignumwpmems_max, type="n",xlim=c(-100,100),ylim=c(0,70),
		xlab="",ylab="",axes=F); axis(1); axis(2,las=2)
		text(jitter(h$s3un100), jitter(h$AS_rignumwpmems_max), h$stateabb,cex=.75,font=2)
		text(jitter(h$s3un100), jitter(h$AS_rignumwpmems_max), h$statenme,cex=.75,font=2)
	mtext("UN Voting Similarity",side=1,line=2.5,font=2,cex=1.25)	
	mtext("Working Party (WP) Member Count",side=3,line=2.5,font=2,cex=1.5)	
	mtext("(Allee Scalera (2011) data)",side=3,line=1,font=2,cex=1)	
	summary(fit <- lm(AS_rignumwpmems_max ~ s3un100, data=h))	
	abline(summary(fit))
	text(90, 0,"(-) 5%", font=2,cex=1.5)
		nobs <-as.data.frame(cbind(h$s3un100,h$AS_rignumwpmems_max))
		nrow(subset(nobs, is.na(nobs[,1])==F & is.na(nobs[,2])==F))
		#45 observations
	mtext("(n=45)",side=3,line=0,font=2,cex=1)
#dev.off()

#pdf("JOP_Appendix_FigureA3_Right_CommitPars.pdf")
	plot(h$s3un100, h$AS_rigcommit_pg, type="n",xlim=c(-100,100),
		ylim=c(0,85),
		xlab="",ylab="",axes=F); axis(1); axis(2,las=2)
		#text(jitter(h$s3un100), jitter(h$AS_rigcommit_pg), h$stateabb,cex=.75,font=2)
		text(jitter(h$s3un100), jitter(h$AS_rigcommit_pg), h$statenme,cex=.75,font=2)
	mtext("UN Voting Similarity",side=1,line=2.5,font=2,cex=1.25)	
	mtext("(Allee Scalera (2011) data)",side=3,line=1,font=2,cex=1)		
	summary(fit <- lm(AS_rigcommit_pg ~ s3un100, data=h))	
		abline(summary(fit))
		text(90, 15,"(-) 10%", font=2,cex=1.5)
	mtext("Commitment Paragraph Count",side=3,line=2.5,font=2,cex=1.5)
		nobs <-as.data.frame(cbind(h$s3un100,h$AS_rigcommit_pg))
		nrow(subset(nobs, is.na(nobs[,1])==F & is.na(nobs[,2])==F))
		#45 observations
	mtext("(n=45)",side=3,line=0,font=2,cex=1)
#dev.off()


#######
# Figure A.4
#	Interaction Models of Geopolitical Alignment and Cold War: 
#	Marginal Effect of Geopolitical Alignment on Time to Apply

####################################################


library(simPH) 
 	# http://www.r-bloggers.com/showing-results-from-cox-proportional-hazard-models-in-r-with-simph/

	#Base sample, time to apply
		aa2 <- subset(d10merge, (ApplyYrI_NotYetApplied_a==1|ApplyYrI_Ind_a==1) 
				&end0>=1 & year>1947	& cty_ind==1); nrow(aa2) #3350

	####################################
	#INTERACTION 1 of 3: COLD WAR X UN VOTING SIMILARITY
	M1 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+s3un100*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		#+old_ptaall_2008extended
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2)
	cox.zph(control.1a)
	cox.zph(M1)
	summary(M1)
	Sim1 <- coxsimInteract(M1, b1 = "s3un100", b2 = "ColdWarPd",
				qi = "Marginal Effect",
                       	X2 = seq(0, 1, by = 1), nsim=1000, spin = TRUE)
		simGG(Sim1, ribbons = TRUE, alpha = 0.5, xlab = "\nCold War",
		      ylab = "Marginal effect of s3un100\n")

		head(Sim1$sims)
	b <- as.data.frame(Sim1$sims)
		head(b)
	cw0 <-subset(b, X2==0)
		nrow(cw0)	#950
	cw1 <-subset(b, X2==1)
		nrow(cw1)	#950
			mean(cw0$QI)
			sd(cw0$QI)
	cw0_cilo_1 <- as.numeric(mean(cw0$QI) - (qnorm(.975)*sd(cw0$QI)))
	cw0_cihi_1 <- as.numeric(mean(cw0$QI) + (qnorm(.975)*sd(cw0$QI)))
	cw0_pe_1	<- as.numeric(mean(cw0$QI))
		cw0_pe_1; cw0_cilo_1; cw0_cihi_1; 

	cw1_cilo_1 <- as.numeric(mean(cw1$QI) - (qnorm(.975)*sd(cw1$QI)))
	cw1_cihi_1 <- as.numeric(mean(cw1$QI) + (qnorm(.975)*sd(cw1$QI)))
	cw1_pe_1	<- as.numeric(mean(cw1$QI))
		cw1_pe_1; cw1_cilo_1; cw1_cihi_1; 

	#INTERACTION 2 of 3: COLD WAR X MEMBER ALLY COUNT
	M2 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+allies_atop_ct_mbrs_zeroed*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		#+old_ptaall_2008extended
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2)
	cox.zph(control.1a)
	cox.zph(M2)
	summary(M2)

	Sim2 <- coxsimInteract(M2, b1 = "allies_atop_ct_mbrs_zeroed", b2 = "ColdWarPd",
				qi = "Marginal Effect",
                       	X2 = seq(0, 1, by = 1), nsim=1000, spin = TRUE)
		simGG(Sim2, ribbons = TRUE, alpha = 0.5, xlab = "\nCold War",
      		ylab = "Marginal effect of allies_atop_ct_mbrs_zeroed\n")

		head(Sim2$sims)
	b <- as.data.frame(Sim2$sims)
		head(b)
	cw0 <-subset(b, X2==0)
		nrow(cw0)	#951
	cw1 <-subset(b, X2==1)
		nrow(cw1)	#950
			mean(cw0$QI)
			sd(cw0$QI)
	cw0_cilo_2 <- as.numeric(mean(cw0$QI) - (qnorm(.975)*sd(cw0$QI)))
	cw0_cihi_2 <- as.numeric(mean(cw0$QI) + (qnorm(.975)*sd(cw0$QI)))
	cw0_pe_2	<- as.numeric(mean(cw0$QI))
		cw0_pe_2; cw0_cilo_2; cw0_cihi_2; # 1.032; 0.9758; 1.088

	cw1_cilo_2 <- as.numeric(mean(cw1$QI) - (qnorm(.975)*sd(cw1$QI)))
	cw1_cihi_2 <- as.numeric(mean(cw1$QI) + (qnorm(.975)*sd(cw1$QI)))
	cw1_pe_2	<- as.numeric(mean(cw1$QI))
		cw1_pe_2; cw1_cilo_2; cw1_cihi_2; # 1.071; 1.016; 1.127

	#INTERACTION 3 of 3: COLD WAR X POLITY
	M3 <-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		#+old_ptaall_2008extended
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2)
	cox.zph(control.1a)
		cox.zph(M3)
		summary(M3)
	Sim3 <- coxsimInteract(M3, b1 = "polity2", b2 = "ColdWarPd",
				qi = "Marginal Effect",
                       	X2 = seq(0, 1, by = 1), nsim=1000, spin = TRUE)

		head(Sim3$sims)
	b <- as.data.frame(Sim3$sims)
		head(b)
	cw0 <-subset(b, X2==0)
		nrow(cw0)	#950
	cw1 <-subset(b, X2==1)
		nrow(cw1)	#950
			mean(cw0$QI)
			sd(cw0$QI)
	cw0_cilo_3 <- as.numeric(mean(cw0$QI) - (qnorm(.975)*sd(cw0$QI)))
	cw0_cihi_3 <- as.numeric(mean(cw0$QI) + (qnorm(.975)*sd(cw0$QI)))
	cw0_pe_3	<- as.numeric(mean(cw0$QI))
		cw0_pe_3; cw0_cilo_3; cw0_cihi_3; 

	cw1_cilo_3 <- as.numeric(mean(cw1$QI) - (qnorm(.975)*sd(cw1$QI)))
	cw1_cihi_3 <- as.numeric(mean(cw1$QI) + (qnorm(.975)*sd(cw1$QI)))
	cw1_pe_3	<- as.numeric(mean(cw1$QI))
		cw1_pe_3; cw1_cilo_3; cw1_cihi_3; 

	min(c(cw0_pe_1, cw0_cilo_1, cw0_cihi_1, cw1_pe_1, cw1_cilo_1, cw1_cihi_1,
		cw0_pe_2, cw0_cilo_2, cw0_cihi_2, cw1_pe_2, cw1_cilo_2, cw1_cihi_2, 
		cw0_pe_3, cw0_cilo_3, cw0_cihi_3, cw1_pe_3, cw1_cilo_3, cw1_cihi_3))	#1.0088
	#min = 0.97
	max(c(cw0_pe_1, cw0_cilo_1, cw0_cihi_1, cw1_pe_1, cw1_cilo_1, cw1_cihi_1,
		cw0_pe_2, cw0_cilo_2, cw0_cihi_2, cw1_pe_2, cw1_cilo_2, cw1_cihi_2, 
		cw0_pe_3, cw0_cilo_3, cw0_cihi_3, cw1_pe_3, cw1_cilo_3, cw1_cihi_3))	#1.0088
	#max = 1.14

#pdf("JOP_Appendix_FigureA4_MargEffects_ColdWar.pdf")
	par(mfrow=c(1,1))
	plot(c(0,1), c(cw1_pe_1,cw0_pe_1), axes=F, type="n", 
	xlab="", ylab="Hazard Ratio",
		ylim=c(.965, 1.20),xlim=c(0,3),
		pch=16, cex=1.5, font.lab=2)
	mtext("Estimated Marginal Effects", cex=1.75, side=3, line=2, font=2)
	axis(2,seq(.95,1.15,.05), las=2,font=2)
		abline(h=1)
	x <- 0.4
	points(x, cw1_pe_3,pch=16,cex=1.25)				#polity
		segments(x, cw1_cilo_3, x, cw1_cihi_3,lwd=2)	#polity
	x <- 0.6
	points(x, cw0_pe_3,pch=16,cex=1.25)				#polity
		segments(x, cw0_cilo_3, x, cw0_cihi_3,lwd=2)	#polity

	x <- 1.4
	points(x, cw1_pe_1,pch=16, cex=1.25)			#unga
		segments(x, cw1_cilo_1, x, cw1_cihi_1,lwd=2)	#unga
	cw1_cilo_1	#0.999
	x <- 1.6
	points(x, cw0_pe_1,pch=16, cex=1.25)			#unga
		segments(x, cw0_cilo_1, x, cw0_cihi_1,lwd=2)	#unga
	x <- 2.4
	points(x, cw1_pe_2,pch=16,cex=1.25)				#ally
		segments(x, cw1_cilo_2, x, cw1_cihi_2,lwd=2)	#ally
	x <- 2.6
	points(x, cw0_pe_2,pch=16,cex=1.25)				#ally
		segments(x, cw0_cilo_2, x, cw0_cihi_2,lwd=2)	#ally

	text(.35, .962, "CW",font=2)
	text(.65, .967, "post-\nCW",font=2)
	text(1.35, .962, "CW",font=2)
	text(1.65, .967, "post-\nCW",font=2)
	text(2.35, .962, "CW",font=2)
	text(2.65, .967, "post-\nCW",font=2)

	text(.5, 1.20, "POLITY", font=4,cex=1.25)
	text(1.5, 1.20, "UN VOTING\nSIMILARITY", font=4, cex=1.25)
	text(2.5, 1.20, "ALLY MEMBER\nCOUNT", font=4, cex=1.25)

	text(.5, 1.175, "Table A.6\nModel (1)", font=2)
	text(1.5, 1.175, "Table A.6\nModel (2)", font=2)
	text(2.5, 1.175, "Table A.6\nModel (3)", font=2)
	
#dev.off()

#######
# Figure A.5
#	Interaction Models of Geopolitical Alignment and Cold War: 
#	First Difference Curves 

######################################################################3

	#PLOT THE FIRST DIFFERENCE CURVE
	#	PERSPECTIVE OF PERCENT FAILURE

	#POLITY ESTIMATES AND GRAPHS

	############
	##TIME TO APPLY FIRST DIFFERENCE CURVE -- POLITY
	##IDENTIFY FIRST DIFFERENCE SAMPLE THAT THE SIMULATION 
	## SHOULD DRAW FROM
	#Model 3
	m3 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		+polity2*ColdWarPd
		+polity2 +s3un100	+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc +gatttotperc
		+lnnew_gdp_67 +lnpcgdp +ColdWarPd 
		+Democratization_MP_ind_Polity + colony_ind_gowa
		+ ptatrade_perc  + PercWorldMbrs + trround_ind 
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2)
		cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE - COMPLETE COUNTRY-YEAR DATA
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1750
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
	datapull<-as.data.frame(consol.1a.sample.in); nrow(datapull)
	affinity	<-as.numeric(summary(datapull$s3un100))
	ally		<-as.numeric(summary(datapull$allies_atop_ct_mbrs_zeroed))
	gatttr	<-as.numeric(summary(datapull$gatttotperc))
	cw		<-as.numeric(summary(datapull$ColdWarPd ))
	gdp		<-as.numeric(summary(datapull$lnnew_gdp_67))
	gdppc		<-as.numeric(summary(datapull$lnpcgdp))
	polity	<-as.numeric(summary(datapull$polity2))
	open		<-as.numeric(summary(datapull$lnopentotperc))
	colony	<-as.numeric(summary(datapull$colony_ind_gowa))
	round		<-as.numeric(summary(datapull$trround_ind ))
	mbrsperc	<-as.numeric(summary(datapull$PercWorldMbrs))
	ptadep	<-as.numeric(summary(datapull$ptatrade_perc))
	democratiz	<-as.numeric(summary(datapull$Democratization_MP_ind_Polity))

	polity.surv0<-data.frame(
		s3un100=rep(affinity[4], 2),	
		allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=c(0,0),			#POST-COLD WAR PD
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=c(polity[2],polity[5]), #move from q1 to q3
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)

	polity.surv1<-data.frame(
		s3un100=rep(affinity[4],2),
		allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=c(1,1),			#COLD WAR PD
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=c(polity[2],polity[5]), #move from q1 to q3
			#polity2=rep(polity[4],2),
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)

#RUN THE SIMULATION 
	store.polity0 <- matrix(rep(NA,68000),nrow=68, ncol=1000)
	store.polity1 <- matrix(rep(NA,68000),nrow=68, ncol=1000)

	for(k in 1:1000){
		datapull<-as.data.frame(consol.1a.sample.in)
		indx<-sample(1:nrow(datapull),5000,replace=TRUE)
		b.data<-consol.1a.sample.in[indx,]
	fitcty	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		#+ColdWarPd*polity2
		+polity2*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		#+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=b.data)
	
		temp <- summary(survfit(fitcty, newdata=polity.surv0))
		d <- data.frame(first = temp$surv[, 1], second = temp$surv[,2], strata = temp$strata, time = temp$time)
		d.subset <- d[d$strata == levels(d$strata)[2], ]
		d.diff <- d.subset[,2] - d.subset[,1]
		store <- rep(NA, 68)
		for(i in 1:68){
			if (i %in% d.subset$time) {
			     store[i] <- d.diff[d.subset$time == i]
			  } else if (i == 1) {
			     store[1] <- 0
			  } else {
			     store[i] <- store[i-1]
			}
			}
		store.polity0[,k]<-store

		temp1 <- summary(survfit(fitcty, newdata=polity.surv1))
		d1 <- data.frame(first = temp1$surv[, 1], second = temp1$surv[,2], strata = temp1$strata, time = temp1$time)
		d.subset1 <- d1[d1$strata == levels(d1$strata)[2], ]
		d.diff1 <- d.subset1[,2] - d.subset1[,1]
		store1 <- rep(NA, 68)
		for(i in 1:68){
			if (i %in% d.subset1$time) {
			     store1[i] <- d.diff1[d.subset1$time == i]
			  } else if (i == 1) {
			     store1[1] <- 0
			  } else {
			     store1[i] <- store1[i-1]
			}
			}
		store.polity1[,k]<-store1
		}
	#started running at:	2:09
	#ended run at:		
	hold0<-store.polity0
	hold1<-store.polity1
	#write.csv(hold0,"FirstDifference,SimulationData,CWintx,Polity,CW0.csv")
	#write.csv(hold1,"FirstDifference,SimulationData,CWintx,Polity,CW1.csv")
		hold0 <- read.csv("FirstDifference,SimulationData,CWintx,Polity,CW0.csv")
		hold1 <- read.csv("FirstDifference,SimulationData,CWintx,Polity,CW1.csv")

#POST-COLD WAR, POLITY INTX
	#TAKE AVERAGES AND 2.5% AND 97.5% OBSERVATIONS
	#CW==0
	ctyest.diff<-ctyci.high<-ctyci.lo<-rep(NA,68)
	for(i in 1:68){
		ctyest.diff[i]<-mean(as.numeric(hold0[i,2:1001]),na.rm=T)
		ctyci.high[i]<-as.numeric(quantile(hold0[i,2:1001],.975,na.rm=T))
		ctyci.lo[i]<-as.numeric(quantile(hold0[i,2:1001],.025,na.rm=T))
		}
	#SAVE THE NUMERIC ESTIMATES TO A FILE
	estimates<-cbind(ctyest.diff,ctyci.lo,ctyci.high)
	estimates0 <- estimates
		#write.csv(estimates0, "FirstDiff,AggregateEst,CWintx,Polity,CW0.csv")
	#USE FIRST LINE IF RAN THE ABOVE; OR, USE .csv FILE WITH SAVED DATA
	#fd0<- as.data.frame(estimates) 	
	fd0<-read.csv("FirstDiff,AggregateEst,CWintx,Polity,CW0.csv")

#pdf("JOP_Appendix_FigureA5_ColdWarIntx0_Polity.pdf")
	plot(1:56, -(fd0$ctyest.diff)[1:56], ylim=c(-.1,.50),xlim=c(0,57),
		ylab="", xlab="",
		type="l", lwd=2,axes=F)
		lines(1:56, -(fd0$ctyci.lo)[1:56], lty=3)
		lines(1:56, -(fd0$ctyci.hi)[1:56], lty=3)
	axis(1,at=c(1,10,20,30,40,50,56),c("1","10","20","30","40","50","56"))
	axis(2,las=2)
	abline(h=0)
	#mtext("POST-COLD WAR", side=3, line=2.5,font=2,cex=1.5)
	mtext("Polity Score", side=3, line=0.5,font=2,cex=2)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Cumulative Probability of Having Applied\nEstimated Difference, 1st to 3rd quartile", side=2, line=3)
#dev.off()


#COLD WAR, POLITY INTX
	#TAKE AVERAGES AND 2.5% AND 97.5% OBSERVATIONS
	#CW==1
	ctyest.diff<-ctyci.high<-ctyci.lo<-rep(NA,68)
	for(i in 1:68){
		ctyest.diff[i]<-mean(as.numeric(hold1[i,2:1001]),na.rm=T)
		ctyci.high[i]<-as.numeric(quantile(hold1[i,2:1001],.975,na.rm=T))
		ctyci.lo[i]<-as.numeric(quantile(hold1[i,2:1001],.025,na.rm=T))
		}
	#SAVE THE NUMERIC ESTIMATES TO A FILE
	estimates<-cbind(ctyest.diff,ctyci.lo,ctyci.high)
	estimates1 <- estimates
		#write.csv(estimates1, "FirstDiff,AggregateEst,CWintx,Polity,CW1.csv")
	#USE FIRST LINE IF RAN THE ABOVE; OR, USE .csv FILE WITH SAVED DATA
	#fd1<- as.data.frame(estimates) 	
	fd1<-read.csv("FirstDiff,AggregateEst,CWintx,Polity,CW1.csv")
		
#pdf("JOP_Appendix_FigureA5_ColdWarIntx1_Polity.pdf")
	par(mfrow=c(1,1))
	plot(1:56, -(fd1$ctyest.diff[1:56]), ylim=c(-.1,.50),xlim=c(0,57),
		ylab="", xlab="",
		type="l", lwd=2,axes=F)
		lines(1:56, -(fd1$ctyci.lo)[1:56], lty=3)
		lines(1:56, -(fd1$ctyci.hi)[1:56], lty=3)
	axis(1,at=c(1,10,20,30,40,50,56),c("1","10","20","30","40","50","56"))
	axis(2,las=2)
	abline(h=0)
	#mtext("COLD WAR", side=3, line=2.5,font=2,cex=1.5)
	mtext("Polity Score", side=3, line=0.5,font=2,cex=2)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Cumulative Probability of Having Applied\nEstimated Difference, 1st to 3rd quartile", side=2, line=3)
#dev.off()

	
#########################
#########################
	
	#ESTIMATE AND PLOT THE FIRST DIFFERENCE CURVE
	#	PERSPECTIVE OF PERCENT FAILURE

	#UN VOTING SIMILARITY ESTIMATES AND GRAPHS

##IDENTIFY FIRST DIFFERENCE SAMPLE THAT THE SIMULATION 
## SHOULD DRAW FROM
	#Model 3
	m3 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		#+ColdWarPd*s3un100
		+s3un100*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		#+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2)
		cox.zph(control.1a)
	summary(control.1a)

		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1750
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014

	datapull<-as.data.frame(consol.1a.sample.in); nrow(datapull)
	affinity	<-as.numeric(summary(datapull$s3un100))
	ally		<-as.numeric(summary(datapull$allies_atop_ct_mbrs_zeroed))
	gatttr	<-as.numeric(summary(datapull$gatttotperc))
	cw		<-as.numeric(summary(datapull$ColdWarPd ))
	gdp		<-as.numeric(summary(datapull$lnnew_gdp_67))
	gdppc		<-as.numeric(summary(datapull$lnpcgdp))
	polity	<-as.numeric(summary(datapull$polity2))
	open		<-as.numeric(summary(datapull$lnopentotperc))
	colony	<-as.numeric(summary(datapull$colony_ind_gowa))
	round		<-as.numeric(summary(datapull$trround_ind ))
	mbrsperc	<-as.numeric(summary(datapull$PercWorldMbrs))
	ptadep	<-as.numeric(summary(datapull$ptatrade_perc))
	democratiz	<-as.numeric(summary(datapull$Democratization_MP_ind_Polity))

	affinity.surv0<-data.frame(
		s3un100=c(affinity[2],affinity[5]),	#move from q1 to q3
			#s3un100=rep(affinity[4],2),
		allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=c(0,0),
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=rep(polity[4],2),
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)

	affinity.surv1<-data.frame(
		s3un100=c(affinity[2],affinity[5]),	#move from q1 to q3
			#s3un100=rep(affinity[4],2),
		allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=c(1,1),
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=rep(polity[4],2),
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)

#RUN THE SIMULATION 
	store.affinity0 <- matrix(rep(NA,68000),nrow=68, ncol=1000)
	store.affinity1 <- matrix(rep(NA,68000),nrow=68, ncol=1000)

	for(k in 1:1000){
		datapull<-as.data.frame(consol.1a.sample.in)
		indx<-sample(1:nrow(datapull),5000,replace=TRUE)
		b.data<-consol.1a.sample.in[indx,]

	fitcty	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		#+ColdWarPd*s3un100
		+s3un100*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		#+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=b.data)
	
		temp <- summary(survfit(fitcty, newdata=affinity.surv0))
		d <- data.frame(first = temp$surv[, 1], second = temp$surv[,2], strata = temp$strata, time = temp$time)
		d.subset <- d[d$strata == levels(d$strata)[2], ]
		d.diff <- d.subset[,2] - d.subset[,1]
		store <- rep(NA, 68)
		for(i in 1:68){
			if (i %in% d.subset$time) {
			     store[i] <- d.diff[d.subset$time == i]
			  } else if (i == 1) {
			     store[1] <- 0
			  } else {
			     store[i] <- store[i-1]
			}
			}
		store.affinity0[,k]<-store

		temp1 <- summary(survfit(fitcty, newdata=affinity.surv1))
		d1 <- data.frame(first = temp1$surv[, 1], second = temp1$surv[,2], strata = temp1$strata, time = temp1$time)
		d.subset1 <- d1[d1$strata == levels(d1$strata)[2], ]
		d.diff1 <- d.subset1[,2] - d.subset1[,1]
		store1 <- rep(NA, 68)
		for(i in 1:68){
			if (i %in% d.subset1$time) {
			     store1[i] <- d.diff1[d.subset1$time == i]
			  } else if (i == 1) {
			     store1[1] <- 0
			  } else {
			     store1[i] <- store1[i-1]
			}
			}
		store.affinity1[,k]<-store1
		}
	hold0<-store.affinity0
	hold1<-store.affinity1
		#write.csv(hold0,"FirstDifference,SimulationData,CWintx,Affinity,CW0.csv")
		#write.csv(hold1,"FirstDifference,SimulationData,CWintx,Affinity,CW1.csv")

	hold0 <- read.csv("FirstDifference,SimulationData,CWintx,Affinity,CW0.csv")
	hold1 <- read.csv("FirstDifference,SimulationData,CWintx,Affinity,CW1.csv")
#TAKE AVERAGES AND 2.5% AND 97.5% OBSERVATIONS
	#CW==0
	ctyest.diff<-ctyci.high<-ctyci.lo<-rep(NA,68)
	for(i in 1:68){
		ctyest.diff[i]<-mean(as.numeric(hold0[i,2:1001]),na.rm=T)
		ctyci.high[i]<-as.numeric(quantile(hold0[i,2:1001],.975,na.rm=T))
		ctyci.lo[i]<-as.numeric(quantile(hold0[i,2:1001],.025,na.rm=T))
		}
	#SAVE THE NUMERIC ESTIMATES TO A FILE
	estimates<-cbind(ctyest.diff,ctyci.lo,ctyci.high)
	estimates0 <- estimates
		#write.csv(estimates0, "FirstDiff,AggregateEst,CWintx,Affinity,CW0.csv")

	#CW==1
	ctyest.diff<-ctyci.high<-ctyci.lo<-rep(NA,68)
	for(i in 1:68){
		ctyest.diff[i]<-mean(as.numeric(hold1[i,2:1001]),na.rm=T)
		ctyci.high[i]<-as.numeric(quantile(hold1[i,2:1001],.975,na.rm=T))
		ctyci.lo[i]<-as.numeric(quantile(hold1[i,2:1001],.025,na.rm=T))
		}
	#SAVE THE NUMERIC ESTIMATES TO A FILE
	estimates<-cbind(ctyest.diff,ctyci.lo,ctyci.high)
	estimates1 <- estimates
		#write.csv(estimates1, "FirstDiff,AggregateEst,CWintx,Affinity,CW1.csv")

#POST COLD WAR PERIOD
#pdf("JOP_Appendix_FigureA5_ColdWarIntx0_Affinity.pdf")
	fd0<-read.csv("FirstDiff,AggregateEst,CWintx,Affinity,CW0.csv")
	plot(1:56, -(fd0$ctyest.diff)[1:56], ylim=c(-.1,.50),xlim=c(0,57),
		ylab="", xlab="",
		type="l", lwd=2,axes=F)
		lines(1:56, -(fd0$ctyci.lo)[1:56], lty=3)
		lines(1:56, -(fd0$ctyci.hi)[1:56], lty=3)
	axis(1,at=c(1,10,20,30,40,50,56),c("1","10","20","30","40","50","56"))
	axis(2,las=2)
	abline(h=0)
	mtext("UN Voting Similarity", side=3, line=0.5,font=2,cex=2)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Cumulative Probability of Having Applied\nEstimated Difference, 1st to 3rd quartile", side=2, line=3)
#dev.off()

		
#COLD WAR PERIOD
#pdf("JOP_Appendix_FigureA5_ColdWarIntx1_Affinity.pdf")
	fd1<-read.csv("FirstDiff,AggregateEst,CWintx,Affinity,CW1.csv")
	par(mfrow=c(1,1))
	plot(1:56, -(fd1$ctyest.diff[1:56]), ylim=c(-.1,.50),xlim=c(0,57),
		ylab="", xlab="",
		type="l", lwd=2,axes=F)
		lines(1:56, -(fd1$ctyci.lo)[1:56], lty=3)
		lines(1:56, -(fd1$ctyci.hi)[1:56], lty=3)
	axis(1,at=c(1,10,20,30,40,50,56),c("1","10","20","30","40","50","56"))
	axis(2,las=2)
	abline(h=0)
	#mtext("COLD WAR", side=3, line=2.5,font=2,cex=1.5)
	mtext("UN Voting Similarity", side=3, line=0.5,font=2,cex=2)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Cumulative Probability of Having Applied\nEstimated Difference, 1st to 3rd quartile", side=2, line=3)
#dev.off()

#########################
#########################
	
	#ESTIMATE AND PLOT THE FIRST DIFFERENCE CURVE
	#	PERSPECTIVE OF PERCENT FAILURE

	#ALLY MEMBER COUNT ESTIMATES AND GRAPHS
	
	##IDENTIFY FIRST DIFFERENCE SAMPLE THAT THE SIMULATION 
	## SHOULD DRAW FROM
	#Model 3
	m3 	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		#+ColdWarPd*s3un100
		+allies_atop_ct_mbrs_zeroed*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		#+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=aa2)
		cox.zph(control.1a)
	summary(control.1a)
		#SAMPLE
			consol.1a.sample <- complete.cases(
				aa2$polity2, aa2$s3un100, aa2$allies_atop_ct_mbrs_zeroed,
				aa2$lnopentotperc, aa2$gatttotperc,
				aa2$lnnew_gdp_67, aa2$lnpcgdp,aa2$ColdWarPd,
				aa2$Democratization_MP_ind_Polity, aa2$colony_ind_gowa,
				aa2$ptatrade_perc, aa2$PercWorldMbrs, aa2$trround_ind)
			sample<-cbind(aa2,consol.1a.sample);consol.1a.sample.in<-subset(sample, consol.1a.sample == TRUE)
		nrow(consol.1a.sample.in)			#obs = 1750
		length(unique(consol.1a.sample.in$ccode))	#countries = 133
		sum(consol.1a.sample.in$ApplyYrI_Ind_a)	#failures = 120
			min(consol.1a.sample.in$year)
			max(consol.1a.sample.in$year) 	#years = 1948-2014
	datapull<-as.data.frame(consol.1a.sample.in); nrow(datapull)
	affinity	<-as.numeric(summary(datapull$s3un100))
	ally		<-as.numeric(summary(datapull$allies_atop_ct_mbrs_zeroed))
	gatttr	<-as.numeric(summary(datapull$gatttotperc))
	cw		<-as.numeric(summary(datapull$ColdWarPd ))
	gdp		<-as.numeric(summary(datapull$lnnew_gdp_67))
	gdppc		<-as.numeric(summary(datapull$lnpcgdp))
	polity	<-as.numeric(summary(datapull$polity2))
	open		<-as.numeric(summary(datapull$lnopentotperc))
	colony	<-as.numeric(summary(datapull$colony_ind_gowa))
	round		<-as.numeric(summary(datapull$trround_ind ))
	mbrsperc	<-as.numeric(summary(datapull$PercWorldMbrs))
	ptadep	<-as.numeric(summary(datapull$ptatrade_perc))
	democratiz	<-as.numeric(summary(datapull$Democratization_MP_ind_Polity))

	ally.surv0<-data.frame(
		s3un100=rep(affinity[4], 2),	
		allies_atop_ct_mbrs_zeroed=c(ally[2],ally[5]), #move from q1 to q3
			#allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=c(0,0),
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=rep(polity[4],2),
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)

	ally.surv1<-data.frame(
		s3un100=rep(affinity[4],2),
		allies_atop_ct_mbrs_zeroed=c(ally[2],ally[5]), #move from q1 to q3
			#allies_atop_ct_mbrs_zeroed=rep(ally[4],2),
		gatttotperc=rep(gatttr[4],2),
		ColdWarPd=c(1,1),
		lnnew_gdp_67=rep(gdp[4],2),
		lnpcgdp=rep(gdppc[4],2),
		polity2=rep(polity[4],2),
		lnopentotperc=rep(open[4],2), 
		colony_ind_gowa=rep(colony[4],2),
		trround_ind=rep(round[4],2),
		PercWorldMbrs=rep(mbrsperc[4],2),
		ptatrade_perc=rep(ptadep[4],2),
		Democratization_MP_ind_Polity=rep(democratiz[4],2)
		)

#RUN THE SIMULATION 
	store.ally0 <- matrix(rep(NA,68000),nrow=68, ncol=1000)
	store.ally1 <- matrix(rep(NA,68000),nrow=68, ncol=1000)

	for(k in 1:1000){
		datapull<-as.data.frame(consol.1a.sample.in)
		indx<-sample(1:nrow(datapull),5000,replace=TRUE)
		b.data<-consol.1a.sample.in[indx,]

	fitcty	<-control.1a<- coxph(Surv(begin0,end0,ApplyYrI_Ind_a) ~ 
		#+ColdWarPd*allies_atop_ct_mbrs_zeroed
		+allies_atop_ct_mbrs_zeroed*ColdWarPd
		+polity2
		+s3un100	
		+allies_atop_ct_mbrs_zeroed
		+ lnopentotperc
		+gatttotperc
		+lnnew_gdp_67
 		+lnpcgdp
		+ColdWarPd 
		+Democratization_MP_ind_Polity
		+ colony_ind_gowa
		+ ptatrade_perc
		+ PercWorldMbrs
		+ trround_ind 
		#+trround_ind:log(end0)
		+cluster(ccode)+strata(Art26_Strata),
		robust=TRUE,data=b.data)
	
		temp <- summary(survfit(fitcty, newdata=ally.surv0))
		d <- data.frame(first = temp$surv[, 1], second = temp$surv[,2], strata = temp$strata, time = temp$time)
		d.subset <- d[d$strata == levels(d$strata)[2], ]
		d.diff <- d.subset[,2] - d.subset[,1]
		store <- rep(NA, 68)
		for(i in 1:68){
			if (i %in% d.subset$time) {
			     store[i] <- d.diff[d.subset$time == i]
			  } else if (i == 1) {
			     store[1] <- 0
			  } else {
			     store[i] <- store[i-1]
			}
			}
		store.ally0[,k]<-store

		temp1 <- summary(survfit(fitcty, newdata=ally.surv1))
		d1 <- data.frame(first = temp1$surv[, 1], second = temp1$surv[,2], strata = temp1$strata, time = temp1$time)
		d.subset1 <- d1[d1$strata == levels(d1$strata)[2], ]
		d.diff1 <- d.subset1[,2] - d.subset1[,1]
		store1 <- rep(NA, 68)
		for(i in 1:68){
			if (i %in% d.subset1$time) {
			     store1[i] <- d.diff1[d.subset1$time == i]
			  } else if (i == 1) {
			     store1[1] <- 0
			  } else {
			     store1[i] <- store1[i-1]
			}
			}
		store.ally1[,k]<-store1
		}
	hold0<-store.ally0
	hold1<-store.ally1
		#write.csv(hold0,"FirstDifference,SimulationData,CWintx,Ally,CW0.csv")
		#write.csv(hold1,"FirstDifference,SimulationData,CWintx,Ally,CW1.csv")
	hold0 <- read.csv("FirstDifference,SimulationData,CWintx,Ally,CW0.csv")
	hold1 <- read.csv("FirstDifference,SimulationData,CWintx,Ally,CW1.csv")

	#TAKE AVERAGES AND 2.5% AND 97.5% OBSERVATIONS
	#CW==0
	ctyest.diff<-ctyci.high<-ctyci.lo<-rep(NA,68)
	for(i in 1:68){
		ctyest.diff[i]<-mean(as.numeric(hold0[i,2:1001]),na.rm=T)
		ctyci.high[i]<-as.numeric(quantile(hold0[i,2:1001],.975,na.rm=T))
		ctyci.lo[i]<-as.numeric(quantile(hold0[i,2:1001],.025,na.rm=T))
		}
	#SAVE THE NUMERIC ESTIMATES TO A FILE
	estimates<-cbind(ctyest.diff,ctyci.lo,ctyci.high)
	estimates0 <- estimates
		#write.csv(estimates0, "FirstDiff,AggregateEst,CWintx,Ally,CW0.csv")
	fd0<-read.csv("FirstDiff,AggregateEst,CWintx,Ally,CW0.csv")

	#CW==1
	ctyest.diff<-ctyci.high<-ctyci.lo<-rep(NA,68)
	for(i in 1:68){
		ctyest.diff[i]<-mean(as.numeric(hold1[i,2:1001]),na.rm=T)
		ctyci.high[i]<-as.numeric(quantile(hold1[i,2:1001],.975,na.rm=T))
		ctyci.lo[i]<-as.numeric(quantile(hold1[i,2:1001],.025,na.rm=T))
		}
	#SAVE THE NUMERIC ESTIMATES TO A FILE
	estimates<-cbind(ctyest.diff,ctyci.lo,ctyci.high)
	estimates1 <- estimates
		#write.csv(estimates1, "FirstDiff,AggregateEst,CWintx,Ally,CW1.csv")
	fd1<-read.csv("FirstDiff,AggregateEst,CWintx,Ally,CW1.csv")

#COLD WAR PERIOD
#pdf("JOP_Appendix_FigureA5_ColdWarIntx1_Ally.pdf")
	fd1<-read.csv("FirstDiff,AggregateEst,CWintx,Ally,CW1.csv")
	par(mfrow=c(1,1))
	plot(1:56, -(fd1$ctyest.diff[1:56]), ylim=c(-.1,.50),xlim=c(0,57),
		ylab="", xlab="",
		type="l", lwd=2,axes=F)
		lines(1:56, -(fd1$ctyci.lo)[1:56], lty=3)
		lines(1:56, -(fd1$ctyci.hi)[1:56], lty=3)
	axis(1,at=c(1,10,20,30,40,50,56),c("1","10","20","30","40","50","56"))
	axis(2,las=2)
	abline(h=0)
	mtext("Ally Member Count", side=3, line=0.5,font=2,cex=2)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Cumulative Probability of Having Applied\nEstimated Difference, 1st to 3rd quartile", side=2, line=3)
#dev.off()

#POST-COLD WAR PERIOD
#pdf("JOP_Appendix_FigureA5_ColdWarIntx0_Ally.pdf")
	fd0<-read.csv("FirstDiff,AggregateEst,CWintx,Ally,CW0.csv")
	plot(1:56, -(fd0$ctyest.diff)[1:56], ylim=c(-.1,.50),xlim=c(0,57),
		ylab="", xlab="",
		type="l", lwd=2,axes=F)
		lines(1:56, -(fd0$ctyci.lo)[1:56], lty=3)
		lines(1:56, -(fd0$ctyci.hi)[1:56], lty=3)
	axis(1,at=c(1,10,20,30,40,50,56),c("1","10","20","30","40","50","56"))
	axis(2,las=2)
	abline(h=0)
	mtext("Ally Member Count", side=3, line=0.5,font=2,cex=2)
	mtext("Years Eligible to Apply", side=1, line=2.5)
	mtext("Cumulative Probability of Having Applied\nEstimated Difference, 1st to 3rd quartile", side=2, line=3)
#dev.off()
##################################################
